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Mathematical Modelling of Compaction and Diagenesis in 


Sedimentary Basins (DPhil Thesis by Xin-She Yang) (Abstract) 


Sedimentary basins form when water-borne sediments in shallow seas are deposited 
over periods of millions of years. Sediments compact under their own weight, causing 
the expulsion of pore water. If this expulsion is sufficiently slow, overpressuring can 
result, a phenomenon which is of concern in oil drilling operations. The competition 
between pore water expulsion and burial is complicated by a variety of factors, which 
include diagenesis (clay dewatering), and different modes (elastic or viscous) of rhe- 
ological deformation via compaction and pressure solution, which may also include 
hysteresis in the constitutive behaviours. This thesis is concerned with models which 
can describe the evolution of porosity and pore pressure in sedimentary basins. 

We begin by analysing the simplest case of poroelastic compaction which in a 1-D 
case results in a nonlinear diffusion equation, controlled principally by a dimensionless 
parameter A, which is the ratio of the hydraulic conductivity to the sedimentation 
rate. We provide analytic and numerical results for both large and small A in Chapter 
3 and Chapter 4. We then put a more realistic rheological relation with hysteresis 
into the model and investigate its effects during loading and unloading in Chapter 5. 
A discontinuous porosity profile may occur if the unloaded system is reloaded. We 
pursue the model further by considering diagenesis as a dehydration model in Chapter 
6, then we extend it to a more realistic dissolution-precipitation reaction-transport 
model in Chapter 7 by including most of the known physics and chemistry derived 
from experimental studies. 

We eventually derive a viscous compaction model for pressure solution in sedi- 
mentary basins in Chapter 8, and show how the model suggests radically different 
behaviours in the distinct limits of slow and fast compaction. When A < 1, com- 
paction is limited to a basal boundary layer. When ’ > 1, compaction occurs 
throughout the basin, and the basic equilibrium solution near the surface is a near 
parabolic profile of porosity. But it is only valid to a finite depth where the perme- 
ability has decreased sufficiently, and a transition occurs, marking a switch from a 


normally pressured environment to one with high pore pressures. 
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Chapter 1 


Introduction 


1.1 Motivation for Modelling Compaction and Diagenesis 


When well-bores are being drilled for oil exploration, drilling mud (a clay suspension 
in water) is used in the hole to maintain its integrity and safety. The mud density 
must be sufficient to prevent collapse of the hole, but not so high that hydrofracturing 
of the surrounding rock occurs. Both these effects depend on the pore fluid pressure 
in the rock, and drilling problems occur in regions where abnormal pore pressure or 
overpressuring occurs, that is in the regions, normally in the sedimentary basins such 
as the North Sea, where pore pressure increases downward faster than hydrostatic 
pressure. Such kind of overpressuring can substantially affect oil-drilling rates and 
even cause serious blowouts during drilling. Therefore, an industrially important ob- 
jective is to predict overpressuring before drilling and to identify its precursors during 
drilling. Another related objective is to predict reservoir quality and hydrocarbon mi- 
gration. An essential step to achieve such objectives is the scientific understanding 
of their mechanisms and the evolutionary history of post-depositional sediments such 
as shales. 

Shales and other fine-grained compressible rocks are considered to be the source 
rocks for much petroleum found in sandstones and carbonates. At deposition, sed- 
iments such as shales and sands typically have porosities of order 0.5 ~ 0.75 or 


50% ~ 75% (Lerche, 1990). When sediments are drilled at a depth, say 5000 m, 
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porosities are typically 0.05 ~ 0.2 (5% ~ 20%). Thus an enormous amount of water 
has escaped from the sediments during their deposition and later evolution. Because 
of the fluid escape, the grain-to-grain contact pressure must increase to support the 
overlying sediment weight. Dynamical fluid escape depends lithologically on the per- 
meability behavior of the evolving sediments. As fluid escape proceeds, porosity 
decreases, so permeability becomes smaller, leading to an ever-increasing delay in 
extracting the residual fluids. The addition of more overburden sediments is then 
compensated for by an increase of excess pressure in the retained fluids. Thus over- 
pressure develops from such a non-equilibrium compaction environment (Audet and 
Fowler, 1992). A rapidly accumulating basin is unable to expel pore fluids sufficiently 
rapidly due to the weight of overburden rock. The development of overpressuring 
retards compaction, resulting in a higher porosity, a higher permeability and a higher 
thermal conductivity than are normal for a given depth, which changes the struc- 
tural and stratigraphic shaping of sedimentary units and provides a potential for 
hydrocarbon migration. 

The compactional fluid escape from the sediments is such a large factor that the 
movement of subsurface fluids must play a dominant role in any attempt to under- 
stand the evolutionary history of geological processes including petroleum formation 
and migration, generation of overpressuring, cementation and dissolution of sedimen- 
tary rocks, fracture formation and dynamical closure, reservoir formation and seals, 
and the formation of ore deposits. Therefore, the determination of the mechanism of 
dynamical evolution of fluid escape and the timing of oil and gas migration out of such 
fine-grained rocks is a major problem. The fundamental understanding of mechani- 
cal and physico-chemical properties of these rocks in the earth’s crust has important 
applications in petrology, sedimentology, soil mechanics, oil and gas engineering and 
other geophysical research areas. 

One purpose of compactional and diagenetic modelling on a basinwide scale is to 
derive an adequate theory to describe the geological processes during compaction, to 
give a series of evolutionary profiles of porosity versus depth, i.e. compaction curves, 


from which geologists and sedimentologists can better understand the burial and 


CHAPTER 1. INTRODUCTION 3 


subsidence histories (Smith 1971). In any attempt to model the dynamic compaction 
of sediments, the main goals are to reproduce, with acceptable agreement and con- 
sistency, four major controls: a) the observed formation thickness, b) the observed 
porosity as a function of depth, c) the observed fluid pressure as a function of depth, 
and d) the observed formation permeability as a function of depth (Lerche 1990). 
These four variables can be calculated in principle from the compaction curves. 

The other important purpose of compactional and diagenetic modelling is to con- 
tribute to a better understanding of how abnormally high fluid pressures come about 
and what factors cause these abnormal pressures to persist for many millions of years 
(Bredehoeft & Hanshaw 1968, Bishop 1979). These high pressures affect seismic in- 
terpretation, mud programs during drilling, and drilling safety. Sediment compaction 
models will be of interest both to the oil industry which always needs better mod- 
els for clay-shale behaviour and to sedimentologists who are concerned with basin 
analysis such as backstripping and burial history. 

The thermal history and the generation of hydrocarbon in a sedimentary basin 
are also closely related to the compaction processes since the thermal conductivity 
and the diagenesis rates depend on the porosity of sediments. The compaction curves 
are also a basis for further studies of petroleum migration. Clay diagenesis is a very 
important process during compacting burial of sediments. Diagenesis is a thermally 
activated reaction in which, for example, water-rich clay mineral smectite dewaters to 
illite, releasing “bound interlayer” water into the fluid system and enhancing the de- 
velopment of overpressuring. Such an illitization process is temperature and pressure 
dependent and is triggered by the catalysis of potassium cations (from K-feldspar). 

In addition the major stage of the smectite-to-illite diagenetic reaction often occurs 
fairly shortly before oz! generation and migration, indicating close organic-inorganic 
interactions. In fact, smectite interlayers may not only incorporate large amount of 
organic products that constitute potential precursors for hydrocarbons, but also act as 
important water reservoirs, that can provide through diagenesis the carrier necessary 
for hydrocarbon migration (Chamley, 1989). Furthermore, it has been recognised 


that overpressuring may often be associated with the formation of seals, which act 
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as barriers to pore fluid expulsion (Hunt, 1990). Within the sealed compartment, 
oil and gas release can build up high pore pressures. The seal formation may in 
turn be related to pressure-enhanced dissolution and reprecipitation of clay minerals. 
All these processes occur in fluid-sediment (water-rock) system, and depend on the 
operating mechanisms of the fluid-sediment interactions. The main purpose of the 
diagenetic modelling is obviously to investigate the operating mechanism of diagenesis 


and reproduce much of the known physics and chemistry of the complex system. 


1.2 Geological Terminology 


The mathematical modelling of compaction and diagenesis is a multi-disciplinary 
study. It is helpful to review the geological terminology related to the present studies. 

Compaction is the process of volume reduction via pore-water expulsion within 
sediments due to the increasing weight of overburden load. The requirement of its 
occurrence is not only the application of an overburden load but also the expulsion 
of pore water. The extent of compaction is strongly influenced by burial history 
and the lithology of sediments. The freshly deposited loosely packed sediments tend 
to evolve, like an open system, towards a closely packed grain framework during 
the initial stages of burial compaction and this is accomplished by the processes of 
grain slippage, rotation, bending and brittle fracturing. Such reorientation processes 
are collectively referred to as mechanical compaction (Kearey & Allen, 1993), which 
generally takes place in the first 1 - 2 km of burial. After this initial porosity loss, 
further porosity reduction is accomplished by the process of chemical compaction such 
as pressure solution at grain contacts. It is worth pointing out that consolidation is 
a term often used in geotechnical engineering and implies the reduction of pore space 
by mechanical loading. 

Diagenesis generally refers to the sum of all those physical, chemical and _bio- 
logical post-depositional modification/reaction processes prior to the onset of meta- 
morphism. Metamorphism is the process of substantial changes to the structure of 


the sedimentary rock by high temperature and pressures. Diagenesis encompasses 
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a broad spectrum of modifications to sediments. Despite its geological importance, 
there is still no universally accepted definition of diagenesis (Rieke & Chilingarian, 
1974). There is no current definitive delimitation of diagenesis either with respect to 
the processes of weathering or metamorphism. In the loosest sense, diagenesis can 
be considered as everything that contributes to making up a sedimentary rock, from 
its weathering near the basin surface to its metamorphism during deep burial. The 
fundamental mechanism behind diagenesis is still less well-documented and there ex- 
ist large discrepancies between laboratory and field data. Diagenesis is influenced by 
burial history, temperature, pressure and pore-fluid chemistry. Diagenesis is dynamic 
as the sedimentary assemblage reacts via the interstitial pore fluids in an attempt to 
equilibrate with the newly established conditions. At diagenetic temperatures and 
pressures, it is very common that the kinetics of diagenetic reactions are slow and 
metastable. Thus, in this sense, diagenesis can be considered simply as low temper- 
ature geochemistry. 

One of the most important diagenetic processes is the smectite-to-illite transforma- 
tion during shale diagenesis. Its reaction mechanism is still under discussion though 
it has received much attention in the last two decades in the literature. One main 
part of our present work is devoted to the mathematical modelling of this important 
diagenetic process. 

Smectite is a family of clay minerals that includes montmorillonite and bentonite 
which is also mainly a kind of montmorillonite-rich clay. The term illite is less a name 
for a definite mineral than a name for a group of substances with composition inter- 
mediate between montmorillonite and muscovite. During diagenesis, montmorillonite 
can release its bounded interlayer-water to form illite which is thermodynamically 
more stable than montmorillonite. 

Dissolution is the diagenetic process by which a solid mineral is dissolved by a pore- 
fluid. There are two fundamental mechanisms for dissolution reactions: transport- 
controlled or surface-controlled dissolutions. The former dissolution reaction is con- 
trolled by the rate of transport of ions to and away from the reacting surface. This 


type of dissolution is typical of fast dissolution by strongly concentrated solutions or 
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of dissolution of highly soluble minerals. In the latter case, dissolution is controlled 
by the reaction rate at the solid-solution interface, and the reaction is relatively slow. 
This type of dissolution is typical of many diagenetic dissolution reactions of relatively 
insoluble minerals in dilute solution with low chemical reactivity. 

Pressure solution/dissolution is the dissolution process under stress. One of its 
most common occurrences is during diagenesis. The increasing vertical load leads to 
dissolution on contact surfaces, and deposition in pore spaces, and thus results in 
(chemical) compaction. The solubility of minerals increases with increasing effective 
normal stress at grain contacts. Pressure dissolution at grain contacts is thus a com- 
pactional response of the sediments during burial in an attempt to increase the grain 
contact area so as to distribute the effective stress over a larger surface. Unfortunately, 
the mechanism and chemistry of the processes are still poorly understood. 

Precipitation is the deposition process of a mineral from a supersaturated pore- 
fluid in either solid form by crystallization or as a gel by flocculation resulting in 
the cementation of the porosity of the host rock. The type of the newly precipitated 
mineral is determined by the type of chemical species in solution and input rate of 
dissolved species into the pore-fluids. Precipitation involves two fundamental pro- 
cesses, nucleation and crystal growth. Nucleation is invariably followed by crystal 
growth, and the two processes are separated by an energy barrier as a result of the 
developing interface between the crystal nuclei and the aqueous solution. Once this 
energy barrier has been surmounted, spontaneous crystal growth, with a net decrease 
in free energy, will proceed until an equilibrium state is achieved when sufficient ma- 
terial is removed from solution so that supersaturation ceases. Experiments show 
that both the rate of nucleation and crystal growth depend upon the supersaturation 
state of the solution. The function of the rate laws is often nonlinear and is poorly 


understood. 
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1.3. Review of Compaction and Diagenesis Models 


Despite the importance of compaction and diagenesis for geological problems, the 
literature of quantitative modelling is not a huge one though the processes received 
much attention in the literature. The effect of gravitational compaction was reviewed 
by Hedberg (1936) who suggested that an interdisciplinary study involving soil me- 
chanics, geochemistry, geophysics and geology is needed for a full understanding of the 
gravitational compaction process. Later in 1959, Weller reviewed the application of 
compaction curves in stratigraphy and structural geology. A more comprehensive and 
detailed review on the subject of compaction of argillaceous sediments was done by 
Rieke & Chilingarian (1974). Audet & Fowler (1992) presented more recently a short 
review on models of compaction. Here we only give a very brief review concerning 
the models of compaction & diagenesis and their developments. 

The mathematical model of compaction and consolidation of shale layers is consid- 
ered as a sediment system consisting of a porous solid phase whose interstitial volume 
is saturated with pore fluid. Due to the action of gravity and the density difference 
between the two phases, the solid phase compacts under its own weight by reducing 
its porosity, thus leading to the expulsion of the pore fluid out of the solid matrix. 
The earliest model about clay consolidation and compaction was proposed by Gib- 
son(1958) based on the earlier work by Terzaghi (1943). This is a linear compaction 
model in which it is assumed that the clay permeability and compressibility are con- 
stant. Gibson’s linear model is sufficiently accurate for modelling thin clay layers 
often encountered in geotechnical engineering. For thick layers and non-constant per- 
meability , the non-linear model was developed by Gibson, England & Hussey (1967) 
and by Gibson, Schiffman & Cargill (1981). Applications of Gibson’s linear model 
investigating sedimentary clay layers were by Bredehoeft & Hanshaw(1968) and Han- 
shaw & Bredehoeft(1968), who modeled diagenesis by considering a layer of source 
rock which produces pore water, leading to overpressuring under the circumstances 
of sufficiently low permeability of the sediments surrounding the source layer. 


The compaction model of shales was developed by Smith (1971) who derived a non- 
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linear compaction model which took into account the dependence of permeability on 
porosity and the dependence of water viscosity on salinity, temperature, and pressure. 
The problem considered was that of a sediment layer growing linearly in time over 
an impermeable basin floor. Water was considered to flow upward or downward 
out of a compacting rock according to Darcy’s law until the pore-water pressure 
within the rock is normal for the depth in question. The conclusion showed that the 
porosity decreases during compaction until a minimum porosity is obtained which 
is determined by the difference between total vertical stress (overburden pressure) 
and pore-water pressure. But Smith’s theory is restrictive in application because the 
compressibility law used by Smith does not include any parameter describing the 
intrinsic strength of the clay sediments. 

The effect and coupling of variation of permeability and temperature with com- 
paction was investigated by Sharp & Domenico (1976) and Sharp (1976) whose results 
are heuristic, but unfortunately technically incorrect. This mistake was finally cor- 
rected by Sharp (1983). Keith & Rimstidt (1985)’s work was similar to the earlier 
work by Smith but the numerical method they used encountered many difficulties 
in the convergence of the numerical results. Therefore, the usefulness of their re- 
sults are restrictive. Bishop (1979) examined a different problem by considering the 
compaction states of thick abnormally pressured shales. The solution predicted the 
interesting characteristic of a density inversion near the overburden shale layer inter- 
face. 

A two-dimensional model was first investigated by Bethke(1985) who investigated 
the case of the temperature dependence of material properties. Unfortunately, as- 
sumptions made in this model are not internally self-consistent, and therefore the 
validity of the results received severe criticism. Some related extensions were anal- 
ysed by Bethke & Corbet(1988) including the porosity dependence of permeability 
and the specific storage. The problem of erosional unloading was treated poroelasti- 
cally by Neuzil & Polluck (1983). An Athy-type constitutive law was used in most of 
these earlier models. 


Audet & Fowler(1992) formulated a rather general mathematical model for the 
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non-equilibrium compaction of clay rocks in sedimentary basins. The model gener- 
alised those of earlier authors. The simplest assumptions were made concerning the 
rheology, but diagensis and thermal coupling were neglected. In this case, their model 
reduced to a generalized consolidation equation, which for the classical Darcy flow is 
a non-linear diffusion equation for the porosity, with a free boundary. The model was 
non-dimensionalized and a robust numerical method was used to solve the non-linear 
diffusion compaction equation. It is interesting that their results only depend on one 
significant dimensionless parameter, the ratio of the Darcy flow rate to the sedimen- 
tation rate. An application of Audet & Fowler’s theory with a detailed parameter 
discussion was made by Audet & McConnell (1992) to investigate the porosity and 
pore pressure evolution for the one-dimensional case in sedimentary basins. Com- 
parison with earlier works shows that the predictions of their model are consistent 
with well data, but it still needs further improvement in the constitutive law for the 
effective stress and for the permeability. 

Wangen (1992) studied the pressure and temperature evolution with a model in 
terms of the void ratio instead of porosity. A new dimensionless parameter is in- 
troduced in this model to characterise the temperature evolution. But the coupling 
between the heat equation and void ratio reduction is a weak one in this model and 
diagenesis is not considered. Luo & Vasseur (1992) investigated the relative impor- 
tance of aquathemal pressuring to geopressure development. This study shows that 
mechanical overloading is the control factor in the development of geopressure but 
the aquathermal effect is less important. Luo & Vasseur’s model and their results are 
similar to Shi & Wang’s model (1986) on pore pressure evolution. Discussion on this 
problem was presented by Miller & Luk (1993) and Luo & Vasseur (1993). 

Diagenesis has been intensively studied in the past two decades, but the attempts 
of dynamical modelling of the process have been made more recently. Several the- 
oretical and computer models have been built. The first model was developed by 
Helgeson (1968) to consider water-rock interactions as a system of coupled dissolution 
and precipitation reactions in which reactions are irreversible and partial equilibrium 


is assumed. Then modified models were proposed by Wolery (1979). These mod- 
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els were then written as software packages PATH (Helgeson, 1968) and EQ3/EQ6 
(Wolery, 1979) and there are later revised versions. The main objections are that 
one has to make an a priori choice of secondary minerals and little information is 
provided on the time-scale of the metastable phases. Other theoretical and com- 
puter models are the REACTRAN model (Ortoleva, Merini, Moore & Chadam 1987) 
and the SOLMINEQ.88 model (Perkins, Kharaka, Gunter & DeBraal 1990). Such 
models predict successive solution compositions and amounts for the dissolved and 
precipitated minerals as water-rock interaction proceeds. As pointed out by Helgeson 
(1979) and Steefel & Cappellen (1990), the assumption of partial equilibrium is only 
justified where the rate of precipitation of a secondary phase is faster then the rate of 
dissolution. However, the precipitation of the stable insoluble minerals may be slow 
even on geological time scales. Therefore, more realistic dissolution and precipitation 
dynamic treatment is essential to diagenetic modelling. Baccar & Fritz (1993) inves- 
tigated a computational geochemical model of sandstone diagenesis and its effect on 
porosity evolution. Their results show that diagenesis effects are very important for 
the evolution of porosity from the point of view of pore fluid chemistry. However, it is 
still difficult to form a clear mathematical model from the existing work on diagenesis. 

Field investigations by Freed & Peacor (1989) in the Gulf Coast and Pearson & 
Small (1988) in the North Sea reveal that diagenesis occurs mainly at burial depths 
from 1 to 2 km in the temperature range from 69° C to 116° C. The illitization of 
smectite with depth in sedimentary basins is observed worldwide and represents one 
of the fundamental reactions in clastic diagenesis. Abercrombie, Hutcheon, Bloch & 
Caritat (1994) analysed the data from oceanic and sedimentary basins and suggested 
that the smectite-illite (S-I) reaction is closely linked to burial parameters such as 
temperature, time and fluid compositions. In a slow sedimentation environment, the 
S-I reaction may begin at temperatures as low as about 50° C, and reach completion 
by about 90° C, while in the rapid sedimentation environment, the S-I reaction may 
not begin at temperatures as high as about 120° C, and reach near completion by 
about 150° C. 


The S-I reaction has received much attention but the nature of both the il- 
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lite/smectite (I/S) mixed-layer and the reaction mechanism are still under discussion, 
and many experiments have been carried out to investigate the kinetic features of the 
S-I conversion (Eberl & Hower 1976; Huang, Longo & Pevear, 1993; Abercromie, 
Hutcheon, Bloch & Caritat, 1994). Thermodynamic analysis shows theoretically that 
quartz and smectite should not coexist at temperatures between 25° and 200°C (Aa- 
gaad & Helgeson, 1983). Lasaga (1984) presented the possible range of activation 
energy variations for a variety of mineral dissolution reactions. Two main mecha- 
nisms have been put forward to explain the S-I reaction process. The transformation 
mechanism suggests that the S-I reaction is a transformation process through mixed- 
layering with (a series of) reordering processes of the intermediate mixed-layer (Hower 
et al, 1976). An alternative modification is a solid-state transformation mechanism 
without mixed-layering. The dissolution-precipitation mechanism involves the pro- 
cesses of smectite dissolution and illite precipitation without mixed-layering. Accord- 
ing to high-resolution electromicroscopic data, the mixed-layering mechanism appears 
to be questionable (Chamley 1989), but Ahn & Peacor (1986) provide a seemingly 
convincing example of a smectite-to-illite transformation rather than a neoformation. 
Although there is no universal consensus, the dissolution-precipitation mechanism is 
theoretically favoured and is consistent with most experimental studies (Chamley, 
1989; Abercromie, Hutcheon, Bloch & Caritat, 1994). 

Potassium cation concentration has an important effect on the reaction rate. At 
is mainly supplied by the dissolution of the K-feldspar. The characterization of K- 
feldspar dissolution rate may be essential for an accurate description of the overall S-I 
process. Four major zones were recognized in the diagenesis of oceanic sediments. The 
apparent lag in illite formation (in deeper zones) may reflect the rate at which the S-I 
reaction proceeds or the availability of potassium for illite formation. The interesting 
correlation with the sedimentation rate and the appreciable variability of the data 
clearly need more systematic work on the mathematical modelling of diagenesis and 
its related consequences. 

The use of correct rate laws are essential to the modelling of water-rock interac- 


tions. Many experiments have been carried out to study the rate laws. However, 
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the laboratory data are not directly applicable to field observations. Unfortunately, 
the discrepancies between field estimates and laboratory measurements of reaction 
rates can be as large as four orders of magnitude (Swoboda-Colberg & Drever, 1993, 
Lasaga et al, 1994). One possibility of explaining this big difference is the effect made 
by coating of mineral surfaces, but this explanation is challenged by the fact that 
extensive etching is widely observed. The fact that diagenesis, which is still imper- 
fectly understood, largely depends on lithology, fluid pressure, geothermal gradient 
and pore fluid compositions is one of the main motivations for us to develop a more 
realistic reaction-transport dissolution-precipitation model in the present work. 

In addition, compaction and diagenesis have been treated separately in conven- 
tional studies. Most available compaction models studied mechanical compaction 
neglecting diagenetic reactions, while the geochemical compaction models mainly in- 
vestigated diagenesis by either prescribing (static) compaction functions or simply 
neglecting the mechanical compaction. 

In summary, the above brief review shows that existing models of compaction and 
diagenesis processes still need more systematical work. In this thesis, we intend to 
extend Audet and Fowler’s (1992) work in the following ways. Firstly, mechanical 
compaction (such as overpressuring), diagenetic (smectite-illite) reactions and ther- 
mal history can be treated simultaneously in a three-dimensional compacting frame 
with a more detailed analysis of some one-dimensional cases of geological importance 
(chapters 2,3 and 4) Secondly, more realistic constitutive models of stress-strain be- 
haviour can be used by employing non-Athy’s type laws and introducing the hysteresis 
during sediment unloading based on experimental data of soils (chapter 5); Thirdly, 
diagenesis can be treated more properly by using more realistic diagenesis models 
such as the first-order dehydration model (chapter 6) and dissolution-precipitation 
model (chapter 7); Fourthly, other related processes such as pressure solution creep 
and fluid geochemistry can also be included in a unified model by utilizing a viscous 
compaction creep law similar to a regelative-flow (chapter 8); Finally, further mod- 
ifications can be developed by considering more realistic basin type and boundary 


conditions (chapter 9). 


Chapter 2 


Mathematical Model 


The general mathematical model of compaction and diagenesis considers the fluid- 
sediment system as a porous medium consisting of multiple mineral species. The 
interstitial volume of the porous solid phase is saturated with pore fluid. Due to the 
action of gravitational overburden load and the density difference between the two 
phases, the solid phase compacts by reducing its porosity, thus leading to the expulsion 
of the pore fluid out of the solid matrix. During compaction and continuous burial, 
the multiple mineral species react and are transported in an evolving pressure and 
temperature environment with a changing rheology. The fundamental underlying 
physical laws to be used are the conservation of mass, the conservation of energy, 
force balance and Darcy’s law. The simple assumptions to be made are related to the 


rheology of the porous medium and the geochemistry of the pore fluid. 


2.1 Audet & Fowler’s Generalised Model for Compaction 


The fundamental model given by Audet and Fowler (1992) can be summarised here 
as follows. Consider a matrix consisting of four interdispersed media: core parti- 
cles (e.g. quartz), two clay minerals, (hydrated) montmorillonite and (dehydrated) 
illite, and free pore water. Let the volume fractions of the respective media (coarse, 


montmorillonite, illite, water) be ¢., dm, ¢;, 1, so that 
Got Omt+ H+ Go =1, (2.1) 
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and suppose that all the solids move with the same averaged velocity u*, while the 
pore water has velocity ul. 


The conservation of mass equations for the four phases will then be in the 


form 
© (pie) + V + (peseu*) = 0, (2.2) 
O 
74 Pmm) an aes (PmPmu’ ) =—Tm; (2.3) 
O 
py Pied +V - (pigiu’) = 7, (2.4) 
«(ob +V- (pda!) = Tw, (2.5) 


where quartz particles are supposed inert, but clay particles can be transformed by 
dehydration processes which release bound water. The rate at which montmorillonite 
is transformed is denoted by 7,,, and this is balanced by a production of illite at rate 
r;, and free pore water at the rate r,,. In fact, the relation r,, = 7; +7T, clearly results 
from the total conservation of mass. 

Diagenesis takes place when montmorillonite ( clay particles with bound water 
between the platelets) releases water to the pore space and is transformed to illite. 
Measured rates of this process in the laboratory (Eberl & Hower 1976) suggest that 
at elevated temperature, this process will proceed very fast from the geological point 
of view. On the other hand, observations suggest that diagenesis is initiated relatively 
suddenly at a temperature 90°C (T.), but then takes place gradually over a depth of 
several hundred meters, which suggests a time scale of the order of a million years. 
This is problematic for the concept of diagenesis as a simple reaction. In fact, the 
mechanism of diagenesis is rather more complicated and is not simply understood. 
Diagenesis may take place via dissolution of montmorillonite in free pore water and 
the subsequent precipitation of silica as illite. Diagenesis is considered here as a one- 
step (first order) dehydration process whose validity is discused in more detail later 
in chapter 7 where a more realistic reaction-transport dissolution-precipitation model 
will be presented. However, the first-order dehydration model is a good approxima- 


tion in the sense of describing the extent of progress of the overall smectite-to-illite 
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transformation without much concern for its detailed geochemical features. Therefore, 


we represent it schematically as 


[clay] - [H2O] (montmorillonite) “5 [clay] (illite) + n|H2O}](free water) (2.6) 


in which we suppose montmorillonite is clay with n moles of bound water per mole 
of clay. From the law of mass action, the prescription of the rates is given by 


nMy 
Mm 


lm = kina Oras y= Ky (57~)PmOm; Tw = k,( 


)PmPm, (2.7) 


where M,,,, M;, My are the respective molecular weights with M,,, = M;+nM,,. The 


reaction rate k, is assumed to follow an Arrhenius law: 


Ea 
k, = Aexp(—Fr)s (2.8) 


where E, is the activation energy which is about 19.6 kcal/mole (Eberl & Hower 1976) 
for the dehydration process, but it may vary in the range of 40-80 kJ/mol (Lasaga, 
1984). R is the gas constant, T is the absolute temperature, and A is a rate factor. 


Let Tp be the surface temperature at the top of the basin; for AT = T — To << To, 


we have 
Eq ioe Te EB, AT 
— Ow (lise), (2.9) 
RT RI Tot AT RH Te 
Hence, k, can also be written as 
Eu Eu 
k= Aexp(— Fr) PS krexp| para (T — To)], (2.10) 
where 
k® = Aexp(— Eo i (2.11) 
, RT 


Denote the heat change per mole during the diagenesis process by AH , and suppose 
that rm, 7; and r,, depend on the temperature and assume that the temperatures of 


each phase are equal, then the energy equation or temperature equation is 


O 
Dp tPeebe oh PHCROR +. PiCiQi ty picidi|T } + ae {[(Pc€ePc TF Prin a picid) Ww 


t+picidu'|T} = V+ (KinVT) — tmAH. (2.12) 


Where c,... are the various specific heats, Ky, is the average thermal conductivity. 


CHAPTER 2. MATHEMATICAL MODEL 16 


According to Fowler (1990), Darcy’s law takes the form 


iG! 3S == (9p da (2.13) 


where j is the unit vector pointing vertically upwards, k is the matrix permeability, 
pl is the liquid viscosity and p! is the pore pressure. 


For a slow flow, the force balance equation can be written 
V-o —pgj =9, (2.14) 


where ois the total stress and the density p = p.@, + pid; (Drew 1983). By employ- 
ing the sign convention for stress in fluid dynamics and using Skempton’s effective 


pressure relation (Skempton 1960, see equation (2.18) in next section) 
—a. =-—o —(l—a)p'd, (2.15) 
the above force balance equation becomes 
V-oe —V{(1—a)p'] — pai = 9, (2.16) 


where o, is the effective stress. 


2.2 Skempton’s Effective Pressure Relation 


Terzaghi (1943) was the first to suggest the principle of effective pressure. According 
to this, the total vertical pressure P at a point in a soil medium consists of two parts. 
One part is carried by water and is continuous and acts with equal intensity in all 
direction. This is the pore water pressure p'. The other part is the pressure carried 
by the soil structure and controls the deformation of the soil structure, and is thus 


called effective pressure pe. Terzaghi formulated this concept as 
pe=P-p, QF) 


which is one of the most important principles of soil mechanics. The modern de- 
velopments on compressibility of soils, shear strength, and lateral earth pressure on 


retaining structures are all based on Terzaghi’s effective pressure concept. Despite 
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its importance, the relation (2.17) is only valid for saturated soils (Skempton, 1960; 
Bear & Bachmat, 1990) . Skempton (1960) extended this relation in a more general 


way, and expressed it in the form 
De = P—(1—a)p’, (2.18) 


where a is a constant. For soils, it may be in the range of 0.1 to 0.5. At pres- 
sures normally encountered in engineering and geological problems, a is very small 
(a <1). Thus, for fully saturated soils, Skempton’s equation degenerates into the 
form of Terzaghi’s equation (2.17) for effective pressure. This corresponds to an 


incompressible and purely cohesive material with a = 0. 


2.3 Constitutive Laws 


2.3.1 Rheological relation for poroelasticity 


The constitutive laws that extend standard linear elasticity to poroelastic materials 
were originally presented by Biot (1941). The constitutive equations were refor- 
mulated by Rice & Cleary (1976) and are most frequently used in the geophysical 
literature. Kumpel (1991) gives a nice review of the poroelastic parameters, and 
more recently Wang (1993) reviews the experimental techniques for measuring the 
static poroelastic moduli and hydrogeologic parameters with particular emphasis on 
the constants that are useful for solving typical geophysical problems. 

Biot’s (1941) linear poroelasticity theory of saturated clay proposed an elastic 
rheological constitutive relation 


2GV 
1— 2p 


Te = 2GeE — ( V-U- p'JI (2.19) 


where € is the strain tensor with e;; = }(0U;/0x; + OU;/0zx;), U is the displacement 
field, I is the second order unit tensor and p, is the effective pressure. G = E/2(1+ 
v),y are shear modulus and Poisson’s ratio respectively, and E’ is Young’s modulus. 


The constitutive relation for porosity is taken to be 


je 7 47-0, (2.20) 
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where @¢ is the initial porosity before deformation. Q and ¥ are two physical con- 
stants. This relation is equivalent to an ordinary elastic medium, with the pressure 


defined by 


1 [BEG =9) 
Seen ea VU: 2.21 
a Q? ~ 30 —2r) een) 
For a saturated clay, Biot suggests Q = oo, y = 1 — a, hence 
3y(1 — 2v) 


d—do =yV-U= (2.22) 


“Ian 
Therefore, pe = pe(¢1) which is an Athy-type law of effective pressure and porosity 
relation. 

For the case of a linear elastic medium, the rheological constitutive relation is 
simplified as 


2 
Oe = 2Ge—(pet+ 3cV -UDI, (2.23) 


with a constitutive relation 


be = —K-V-u’, (2.24) 


where’ denotes d/dt, = a +u*-V and K, is a constant. To follow o, with a material 


element, we have 
doe 


dt, 


It is worth pointing out that the rheological equation of state should be objective. 


= 2Gé — (p. + “Gv uw )L. (2.25) 


That is to say, the rheological relation of stress-strain should be invariant under the 
coordinate transformation. This is not always guaranteed due to the complexity of 
the rheological relations (Bird, Armstrong & Hassager 1977). Fortunately, for one- 
dimensional zrrotational flow, the equation is invariant and all the different equations 
in corotional and codeformational frames degenerate into the same form. In the 
one-dimensional case we will discuss below, we can take this for granted. 

For the very simple case of a one-dimensional model, the effective stress tensor is 
in the form of 


Ce = diag(—o1, —01, —03). (2.26) 


By using equation (2.23), we have 


)De- (2.27) 
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2.3.2 Permeability 


The permeability, for the convenience of later usage, can be expressed in the normal- 
ized form (Smith 1971) 


ee Cen eos (2.28) 


o 


where ko and @p are the permeability and porosity at the top of the basin. m is a 
positive number which characterizes how quickly the permeability decreases as the 
porosity is reduced. The unit of permeability is Darcy (1D = 10~'’m?). Typical 
values of the permeability of clays are in the range of 1.5 x 107-8 ~ 1.5 x 1073 Darcy 


for the porosity range from 0.33 to 0.8 (Lambe & Whitman 1979). 


2.3.3 Thermal conductivity 


Thermal conductivities of sedimentary rocks vary with porosity. High-porosity uncon- 
solidated rocks have low values of thermal conductivities, while nearly fully compacted 
sediments with low porosity have high values. To calculate the averaged thermal con- 
ductivity Ky, of a porous medium, we use a rough quasi-empirical relation (Lewis & 
Rose, 1970) 

Kin = Kol), (2.29) 


s 


where K, is the thermal conductivity of pore water and K, the thermal conductivity 
of sediment matrix. Ko and ¢p are the thermal conductivity and porosity at the top 


of the basin. 


2.4 One-dimensional Model 


In order to simplify the following calculations and to compare the results with earlier 
work, we will consider a one-dimensional compaction model in a basin b(t) < z < h(t) 
(Fig. 2.1), where h is the ocean floor and b is the basement rock, instead of more 
general cases in two or three dimensions. This one-dimensional compaction model is 
applicable to the case in which the basin depth is small compared to its length and 
width. 


CHAPTER 2. MATHEMATICAL MODEL 20 


| im, 


Z Ocean floor: z = h(t) 


0 Basin basement: z = b(t) 


SSS RAG ASA ARAASS 


Figure 2.1 One-dimensional compacting sedimentary basin. The coor- 


dinate z is directed upwards. 


2.4.1 1-D governing equations 


For convenience in the following discussion, we put p! = p. We will investigate the 


simplest behaviour of non-linear compaction restricting our attention to the case 


where the solid species have density pp = Pm = Pj = Ps =constant and specific heat 


Co = Cm = CG = Cs =constant. With these simplifications, we can easily obtain the 


governing equations from the above section. 


Mass conservation 


ots += 5 Pou" Ve (2.30) 
9m O 
Oo ee BAM 
Men 2 ut) = (eT bn (2.33) 
bdeotomt hth =, (2.34) 
Darcy’s law 
k Op 


gi(u! — u®) = —-(= + pig), (2.35) 
Ll 
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Energy conservation 


© (o.ca(l — $1) + pergilT} 


O O OT 
+7 i lPscs(1 — pur + picipiu |T} = 9; King) — kppsbmAd, (2.36) 
Force balance 
O 4G 
De! (14 3K pe — (1 — a)p] — [ps(1 — or) + pidilg = 9, (2.37) 
Constitutive relation 
Pe = De(d1)- (2.38) 


These are nine equations for nine unknown variables: four for volume fractions 


be, Pm, i, b1, two for velocities u*, u!, 


one for temperature 7’, and two for effective 
pressure p, and pore water pressure p. 
In order to get an expression for u*, we add the four equations of mass conservation 


together and thus have 


M, 
Z “\obm with 6 = ps/ pr. (2.39) 


O i Ss) __ 
7, [ere + (1 — j)u*] = kr(d — 1)( M,, 


By using Darcy’s law, the above equation becomes 


OW. Ok Op 


a =f ae 2.4 
Fe = Belge + Pua + Bl EI (2.40) 
Integrating z from 0 to z, we obtain 
Op, 
v= nee + pig) + (6 — 1)( “krmdz + b(t), (2.41) 
and 
dul + (1 — due = (6 — 1)( * ke@mdz + b(t). (2.42) 


2.4.2 Boundary conditions 


The related boundary conditions for the nine governing equations are as follows. If 
we take b(t) as a known boundary, but h(t) as unknown, then we still require bound- 
ary conditions on wu!,u°,p,p.,T' for the equations. Obviously, the natural boundary 
conditions are the following: 

boundary conditions at z = b: 


ud =u! = 6; (2.43) 
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a kinematic condition at z = h: 
h=m,.+u', (2.44) 
where m, is the sedimentation rate at z = h. Also at z= h, 


bi = bo = do, (i-e.,Pe = 0), P= Do (2.45) 


and 
Ge = G0, Pi = Pin, Om = Pmo, (2.46) 


where po is the overburden pressure, e.g. due to ocean depth. ¢.09, ¢i0, 6mo and ¢o 
are the values at the top of basin during sedimentation. 


The boundary conditions for the temperature (or energy) equation become 


Tt=0,2=h)=T7> and a(h2 = 8) = — (2.47) 
z 


where qo is the heat flux at the bottom of the basin. This corresponds to a constant 
temperature 7p at the top of the basin and a constant heat flux at the base. Equation 
(2.44) gives the moving boundary h(t), and therefore we have the number of conditions 


which the equations require. 


2.5 Non-dimensionalization 


We define a length-scale d by writing 


Pe = (Ps — pr)gdp( Gu), (2.48) 


and require that p = O(1). Meanwhile, we scale z with d, u® with m,, time t with 
d/mm;, pore pressure p with (p,— 1)gd, permeability k with ko, heat conductivity Ky, 
with Ko, temperature T with qod/Ko, k, with k° and AH with qo/(msps); thus we 
have 

k = kok, kp = kokp Kin, = KoK, (2.49) 


and 


d 
Parnes, a 


AH. 2.50 
Ko MsPs ( ) 
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The dimensionless form of equations (2.30)-(2.37) is then 


Mass conservation 


a _ 
ba a < (on! = Rk,a10Om, 


e+ Omt H+ = 1, 


Darcy’s law 


Energy conservation 


O 
DE la(1 — ¢) + dO} 


0 0, ~00 


tay tlatl — bial + yu'|O} = i ae ARK pdm AH, 


Force balance 


where 


and define 
y= | krdmdz, clearly i =0 on z=0. 
0 
For the diagenesis parameter, we have 


Eaqod 


Kp = exp((O) with B = RKoTe’ 


23 


(2.51) 
(2.52) 
(2.53) 


(2.54) 


(2.55) 


(2.56) 


(2.57) 


(2.58) 


(2.59) 
(2.60) 


(2.61) 


(2.62) 


(2.63) 


where O = (T — 75) Ko/qod is the dimensionless temperature with reference to the 


surface temperature To. 
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The dimensionless expressions for u* and u! from (2.41) and (2.42) now become 


ue = a(S +r) +(6—1)6Ry +), (2.64) 
du! + (1 — d))u® = (6 — l)ayRwy + 6. (2.65) 


The boundary conditions in the dimensionless form are 


Bo) 2 —(1+ar—¢;)=0at z=6), (2.66) 
(l1—a)h=(1—a)m— Mile (o) Se + (1+ar— ¢,)| 
+(1-a)(6- (A )RY + (1 -apb Pi (2.67) 
d= 0, be = Geo, Pi = Pi0, Pm =Pmo at z=h, (2.68) 
666s So aad (tz =) =-=. (2.69) 


Here, 7 is the dimensionless sedimentation rate which is 1 if it is constant, or O(1) 
if time-varying. 

It is very interesting that the above derived dimensionless porosity, temperature 
and diagenesis equations are based on eight dimensionless parameters. The five pa- 
rameters r, a, 6, a1, @ are constants to some extent. The other three parameters, 
namely A, A, R, are the governing parameters controlling the whole evolution pro- 
cess. It is worth pointing out that the parameters r and 6 are not independent. 
A = ko(ps — pdg/ums, A = Ko/picmsd and R = k°d/m, are parameters which char- 
acterize the porosity, the temperature and diagenesis evolution, respectively. Here ko 
is the permeability at the top of the basin, jz is the viscosity at the top of the basin, 
pi, c, are density and the heat capacity of fluid (water). The parameter Ko in A is 
the bulk heat conductivity of the sediments at the top of the basin. 


2.6 Determination of Model Parameters 


It is useful for the understanding of the solutions to get an estimate for A, A and R by 
using values taken from observations. The model parameters are chosen by referring 


the values given by other authors (Smith 1971, Sharp 1976, Sharp & Domenico 1976, 
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Eber! & Hower 1976, Bethke & Corbet 1988, Lerche 1990, Audet & Fowler 1992). The 
values used in the present model are d ~ 1 km, ky ~ 1 x 10°’ m?, p, ~ 2.6 x 10° ke 
m3,g~10ms-?, p, ~ 1x10? kg m3, pp ~1x10-3Nsm?, m, ~ 300m Ma! = 1x 
10-"' m s-!, e, ~ 500 J Kg? K-1,¢, ~ 4200 J Kg"! K-!, Kp ~1x15Wm! K-}, 
Ty ~ 280 K, T, ~ 363 K, Ey ~ 8.18 x 104 Jmol“! and k® ~ 1 x 107s"! ; then 
Bx 2.3,’\ 1, AX 30 and R + 0.01. Therefore, \ = 1 defines a transition between 
the fast sedimentation (A << 1) and slow sedimentation (A >> 1). The parameter A 
, which is the ratio between the permeability and the sedimentation rate, governs the 
evolution of the pore pressure and porosity in sedimentary basins. High sedimentation 
rate may gives rise to excess pressures even in the basins with moderate permeability. 

Similarly, the parameter A also defines a transition. A << 1 shows that the 
temperature solution is dominated by the constant growth of the basin thickness due 
to fast sedimentation, while A >> 1 shows that the sedimentation rate has little 
influence on the temperature solution. The parameter R characterizes the effect of 
diagenesis on compaction. 

An initial porosity of ¢@) = 0.5 for pore water at the top of the basin has been used 
by other authors (Smith 1971, Sharp 1976, Bethke & Corbet 1988, Audet & Fowler 
1992). Initial porosity 0.2 for montmorillonite, 0 for illite and 0.3 for quartz are used 


in the following computations. 


2.7 Overpressure Definition 


The hydrostatic pressure at z is defined as 


A(t) 
Ph = | pigdz. (2.70) 


The overburden pressure at z is defined as 


h(t) 
P= [ (= i)ps + drpilgde. (2.71) 


The excess pore pressure or abnormal overpressure p, is defined as 


Pa = P— Ph; (2.72) 
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which is the pressure in excess of the hydrostatic pressure. 
By using these definitions and employing the force balance equation ( 2.58), the 


dimensionless differential forms of the above definitions are 


OF 
Oph = 
ODa a Op 
(1 —a) ieee (1+ ar — ¢)). (2.75) 


It can be seen that pressure profiles can be easily calculated from the porosity profile 
or compaction curve of @ versus z. Therefore, the main target is to find evolving 


features of the compaction curves. 


Chapter 3 


Numerical Simulations 


3.1 <A Simple Case 


In a moving frame of reference, it is obvious that b = 0 can be selected, equations 
(2.52), (2.53), (2.54),(2.56), (2.57) and (2.58) then form a free boundary problem 
for dm, o;, ¢; and ©, depending essentially on three parameters A, A and R. For 
simplicity we also take AH = 0 in these equations. Based on the work of Smith 


(1971), Sharp (1976) and Audet & Fowler (1992), we adopt the following constitutive 


functions: 
P = In(¢o/¢1) — (¢0 — 61); (3.1) 
k = (1/40), m=8, (3.2) 
K =(K,/K,)*"™, K,/K, = 0.3, (3.3) 
m=1, k, =exp(@@). (3.4) 


By using these constitutive relations together with the force balance equation, and 
eliminating u®, u! and p in equations (2.56), (2.58) and (2.64), we can the obtain 
coupled non-linear diffusion equations for ¢.,@m, @; and © whose forms are suitable 
for numerical calculations and asymptotic analysis. These equations are 


Equations for volume fractions 


Ob 


: F) 
(1-0) = 0 fio — dnl et ~ (14 


Bon eg 


27 
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~(1 ~a)(8— 1a RULE + Foe), (3.5) 
Obm _ = O + 1 0¢, _ ar 
(1 — aa x = = a)RKPm — Ag, tkeom(d _ ONE Oz = (1 I fuck ral 
—(1 ~ (6 — 1) RW LO” + Fg), (3.6) 
Oh — , Ope 9,1 Ob; ar 
(=) = AS (ha - oS - 04+ yy 
+(1—a)Rhpaidm ~(1—a)(6— 1arRWSE + (1 )omh (37) 
Temperature equation 
ele) 0,00 2 
(1 — a)[a(1 — ¢)) + re =(1- aJAn (Ka) — (1—a)Rk,a,(6 — a)bmO 
00 ~ 1 0¢, _ ar ,,00 
—(6-1)(1- a)arba— —(a—1)Ak(1 -— 15 geet (14 = 7a Dy? (3.8) 
where 
as ko(ps = pig _ Pt . (3.9) 
UMs Ps — Pl 
Bo AIOE, pesca (3.10) 
picim,d pic 
_ kd mMy ._ ps 
= ze ae oars (3.11) 
and 
The related boundary conditions (2.66)-(2.69) become 
Od] ard, = 00 =) _ ol _ 
eg ae ee ie 
and 
e 10 
8,(t,h(t)) = do, (1— adh = (1 arin + ARCA ~ OF" ~ (14+ 
+(1—a)(d — 1l)a,yRy at z =h. (3.14) 


where 6; = $c, i,m, or oy; and ¢. + ¢; + dm+¢, = 1. The condition ¢)(t) = ¢o at 
z = h(t) is equivalent to p(t) = 0 (the effective pressrue is zero). 


The dimensionless form for the excess pressure (2.75) is then 


(1 — a) oP = (1 - ¢)(—=—- - 1) -ar, (3.15) 
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with a bounday condition p, = 0 at z = h(t). 

It is based on these equations that the moving boundary problem will be solved 
numerically by using the predictor/corrector implicit finite-difference method pre- 
sented by Meek & Norbury (1982), which is very robust for the non-linear parabolic 


equations. 


3.2 Audet & Fowler’s Case 


If we set a = 0, R = 0 (no diagenesis) and leave out the temperature equation in 
the previous section, we then get a very special case which was considered by Audet 
& Fowler (1992). The equation for ¢; degenerates simply to a general non-linear 


diffusion equation 


“en =15 Rd = 61)*15 - a} (3.16) 
h= 1+ ARC = 65S — 1), (3.17) 
with boundary conditions 
a me eT z=0, (3.18) 
Oz 
Olt.) = o¢al 2= nh. (3.19) 


which was discussed in detail by Audet & Fowler (1992). 


3.3. Finite Difference Approach 


In order to solve the highly coupled non-linear equations in this work, an implicit 
numerical difference method is used (Smith 1985). The essential equations describing 
for porosity and temperature are of the standard non-linear parabolic form (Meek & 
Norbury, 1982) 

ty =F (GEA tie F(t) (3.20) 


n+1/2 


The first stage gives u as a solution of the following equation 


ee ‘ 
A mp2 = ur) = ( )F (xj, 0°78? ur) e2urt/? 


Az? 
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L 
cm) ite Ap aA 
ey Uh 


where 62u; = (uiy1 — 2u; + ujz_1) and 6,u; = (1/2)(ui41 — uj_1). The second stage 


bx U; ), (3.21) 


gives u?*? as a solution of the following equation 
1 nt+1 n l n+1/2 , nt+1/2) 27. n+1 n 
Ay i — uz) = (Sagar ut » Uy )oz(ug™ + uz) 


1 


+f (xj, t?tV?, Ayo 2); (3.22) 


The convergence is second-order in space for this method, and O(At)?~< in time, 
where ¢€ is a small number less than 1/2. 

The computational convergence of the calculation of this method has been tested 
by 1) changing the grid number from 5 to 1000 in space and from 10 to 5000 in 
time, and by 2) comparing with the results of asymptotic results. The changes of 
grid intervals all result in the same converged results which conform well to the 
asymptotic solutions. This shows that this method is robust for the solution of the 


equations encountered in our problems. 


3.4 Numerical Results 


By using the above mentioned implicit numerical method, we can solve the equa- 
tions numerically for various values of A, A, R and @. We used a normalized grid 
parameterized by the fixed domain variable Z = z/h(t). This will make it easy to 
compare the results of different times and different depths with different values of 
dimensionless parameters in a fixed frame. This transformation maps the basement 
of the basin to Z = 0 and the basin top to Z = 1. The numerical method was first 
tested in MATLAB and later transformed to FORTRAN codes with double preci- 
sion. The calculations were mainly implemented for the time evolutions in the range 
of t = 0.5 ~ 10 since the thickness in the range of 0.5km ~ 10km is the one of interest 
in the petroleum industry and in civil engineering. Preliminary numerical results are 


presented and explained briefly below. 
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3.4.1 Comparison with Audet & Fowler’s results 


In Audet & Fowler’s case (a = 0 and R = 0), we get a moving boounadry problem 
(3.16) (Section 3.2). Solving this problem with different values of \ and time t, we 
have the following numerical results. 

Porosity evolution with different \ values, corresponding to different sedimen- 
tation rates, are calculated. The numerical results are shown in Fig. 3.1 for different 


values of A = 0.01, 0.1, 1, 10, 100, at a fixed time t = 5. 


Porosity Evolution and Transition in Compaction 
T T T T 


0.97 t=5, a=0 


0 01 02 ak 03 
Figure 3.1 Porosity evolution with different values of \ or sedimentation 
rates. Z is scaled height, and the different values of \ are given along 
the curves. This figure shows that porosity evolution is essentially con- 
trolled by . A porosity boundary layer develops near the basement in 
a rapid sedimentation environment (A = 0.01) while porosity decreases 
nearly exponentially in a slow sedimentation environment (A = 100). 

For the case of A = 100 in Fig. 3.2 and A = 0.01 in Fig. 3.3, different results 
for different evolution times are plotted as porosity versus depth. Fig. 3.2 and Fig. 
3.3 are essentially the same results as those discussed by Audet & Fowler (1992). It 
is clearly seen in Fig. 3.2 that there exists a travelling wave solution of porosity @ 
for the large \ case in the top region where the porosity profile is only a function of 
depth z — h(t). On the other hand, a boundary layer develops near the basement in 
Fig. 3.3 for the case of small i. 
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Porosity versus depth at various times 


z-h(t) 
w 


0 0.1 0.2 0.3 0.4 0.5 
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Figure 3.2 Porosity versus depths at various times with a fixed value of 


A = 100. z — h(t) is the depth measured from the basin top. 


Porosity versus depth at various times 
T T 


0 
-0.5+ a=0, R=0, lambda=0.01 | 
tel 
-1b 4 
1st 
t=2 
-2b 
E25 
& 
t=3 
3b 
3.5¢ 
t=4 
-4b 
4.5 
t=5 
5 
0.4 0.42 0.44 0.46 0.48 0.5 


Porosity 


Figure 3.3 Porosity versus depths at various times ¢ for \ = 0.01. z— 


h(t) is the depth measured from the basin top. 
The results in Fig. 3.1-3.3 show that the parameter is the most important 


dimensionless parameter controlling the degree of compaction and overpressure. In 
the case of high permeability and low sedimentation rate (A >> 1), the pore water 
will leave the sediments at almost the same rate as the increase in the overburden 
load. The sediment column will remain nearly hydrostatic and the compaction will 


be almost maximal. Thus, porosity decreases nearly exponentially in the top region. 


CHAPTER 3. NUMERICAL SIMULATIONS 33 


While, in the opposite case, with low permeability and high sedimentation rate (A << 
1), the water is nearly unable to escape from the sediments at the same rate as the 
increase in the overburden load. Water gets trapped in the pores, water pressure 
builds up, and the compaction is very small. This results in a porosity boundary 
layer near the basement. 

Effects of a on porosity are computed. Fig. 3.4 gives the results of different 
values of a = 0,0.3,0.6,0.9 for the same evolution time t = 5 with values of \ = 
1, A = 1 and R = 0. This clearly shows that a has a significant effect on porosity. 
In the extreme case, a = 1, which corresponds to the elastic perfectly-compacted rock 


sediments, the porosity will be zero at all depths. 


Effect of a on Porosity 
T T 


0.97 t=5 


0 0.1 0.2 0.3 0.4 0.5 
Porosity 


Figure 3.4 Effect of a on porosity for \ = 1 and t = 5. Z is scaled 
height, and the different values of a are given along the curves. This 


clearly shows that a has a significant effect on porosity evolution. 


Basin thickness is calculated for different values of 4 = 0.1,1,10 for the case of 
t < 10 with all the other values fixed as before. 

The results in Fig 3.5 demonstrate that the thickness of the basin for different 
sedimentation rates is always nearly linear, but the slopes can vary (0.98, 0.68, 0.59 


respectively for A=0.01, 1, 100). 
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Figure 3.5 Thickness h(t) versus time t at different sedimentation rates 


(A = 0.01, 1, 100), with a initial value of h(0) = 0. h(t) increases 


Thickness with time for different sedimentation rates 
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Subsidence and fluid flow velocity are presented in Fig. 3.6. The solid velocity 


at the top of the basin shows the compaction-driven subsidence velocity of the basin 


top. The fluid velocity u! and Darcy velocity (dashed curve) shows the compaction- 


driven fluid flow fields at different depths. 


Figure 3.6 Subsidence velocity u® (solid), fluid velocity u! (solid) and 
Darcy velocity (u! — u’) (dashed) versus scaled height Z at time t = 5 


for \= 1. 


Fluid flow and Darcy velocity 
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3.4.2 The development of excess pressure 


The overburden (P), hydrostatic (p;,) and pore (p) pressures are calculated for two 
cases of A = 0.01 and AX = 100. The results are shown in Fig. 3.7 and Fig. 3.8, 
respectively. 

Figure 3.7 shows that the water is almost unable to escape from the sediments at 
the same rate as the increase in the overburden load in the case of low permeability 
or high sedimentation rate (A << 1). Fluid gets trapped in the pores, pore water 
pressure builds up, and the compaction is very small. The excess pressure develops 
proportionally to basin thickness. 

Figure 3.8 shows that pore water will leave the sediments at almost the same 
rate as the increase in the overburden load in the case of high permeability or low 
sedimentation rate (A >> 1). The sediment column will remain nearly hydrostatic 
and the compaction will be almost maximal. Excess pressure does not occur for short 


times or in the top region but may develop at large times in the lower region. 


scaled height: Z 
j=5 
wn 


0.47 overburden 


0.37 hydrostatic 


i pore 
0 0.5 1 1.5 2 
scaled pressure: p/h(t) 


Figure 3.7 Hydrostatic, pore and overburden pressures at t = 5 for the 
case of \ = 0.01. Because water is almost unable to escape from the 
sediments at the same rate as the burial, water gets trapped in the 


pores and pore pressure builds up. 
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scaled height: Z 
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Figure 3.8 Hydrostatic, pore and overburden pressures at t = 5 for the 
case of X = 100. Excess pressure does not occur for short times or in 


the top region but may develop at large times in the lower region. 


3.4.3 Temperature evolution 


For simplicity, we now put R = 0, a = 0, @ = @, then the governing equation for 


temperature evolution becomes 


OO 0,00 
is Se Be! he 
(a — )AK(L — "(SF — SE (3.23) 
with boundary conditions 
OO 1 
and 
Q=0 ee a eee eee =h (3.25) 
—_— “0; = @ Oz at z= fh. : 


This is a moving boundary problem which can be solved numerically. 


Temperature profile for values A = 0.1,1,10 is shown in Fig. 3.9 with other 


parameters fixed (t=5,a=0,a=0.3,A=1, R=0). 
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Temperature Evolution & Comparison with Analytical Solution 
T T T T 
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Figure 3.9 Comparison of analytical solutions (dashed) for temper- 
ature evolution with numerical results (solid) for different values of 
A =0.1, 1, 10. Z is scaled height. The asterisk (*) corresponds to the 


steady-state solution. 


The non-linearity seen in Fig. 3.9 for A = 0.1 is due to the effect of the fast mov- 
ing boundary and the conductivity function K (¢,). The calculations show that the 
heat convection term in equation (3.23) has no significant effect on the temperature 
evolution. This result is in accordance with other authors (Bethke, 1985, Deming, 
Nunn & Evans, 1990) who have pointed out that convective heat transfer is less im- 
portant in the one-dimensional compaction models, but may be important in two- or 


three-dimensional models with lateral fluid flow. 


3.4.4 Heat conduction with constantly moving boundary 


From the numerical simulations, we see that the heat convection term has no signif- 
icant effect on the changes of temperature evolutions. This can be understood from 
the fact that the second term on right hand side of equation (3.23) is equivalent to 
(a—1)(1—¢)u'00/dz. This means that there is no significant difference in the tem- 
perature profiles when a changes from 0.3 to 1. The analysis in the next chapter will 
show that the convective term is O(A) << 1 for slow compaction (A << 1), while for 


fast compaction (A >> 1), Athy’s solution suggests that ¢,/¢ © 1, then the convec- 
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tive term is also small compared to the conduction term. In addition, the numerical 
results in Section 3.4.1 show that h(t) almost linearly depends time t. Therefore, we 
will mainly concentrate on the solutions of the standard equation of heat conduction 
with a constantly moving boundary. To approximate this, let h(t) = Ut. For \ << 1, 
U = 1; for \ >> 1 then U & 0.59. To simplify the analysis, we assume a = 1 and 


K =constant, then the temperature equation becomes 


00 OO 
a = AK Se (3.26) 
with 
Olt. z= h):=—0-and tt e=0) ==. (3.27) 


This is a diffusion equation with a specified moving boundary. The solution of this 
problem can be constructed by employing Green’s function method and using the 
method similar to Gibson’s (1958) approach for a consolidation problem with a con- 


stantly increasing thickness. We assume the solution of the following form 


Oa. t= 2(AKi)' ierfc| é 
kK 2(AKt)!/2 
A) 00 =? 2 
Sar eV ee g(Cexe(- = oa (3.28) 
where ierfc(¢) = eee — Cerfe(¢), and erfc(¢) = 1— = Se Ig e =n? 


It is easy to check that the above solution satisfies the temperature equation and 
the boundary condition at z = 0. Therefore, we are at liberty to regard g(¢) as an 
arbitrary function which must be chosen to satisfy the upper boundary condition at 
z= h(t). This requirement is met if g(¢) satisfies 

h2 


(4r)'/? At ierfe[ U, Jexp(——) 
a AAKt 


2-2 
=f 9 sinh sl, we exp(— TS Ja. (3.29) 


By substituting h(t) = Ut, changing e eae. Cah t= 1/4A Kr and using the 
Laplace integral technique, we obtain 


By usuerre ee 
g(¢) = on ee 
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1 | i fel( Oi U? 
Sees Sehr 16A2r7 


Substituting this back into the solution, we will have an integral form of the solu- 


)/Jexp/ + 7¢*)\dr. (3.30) 
tion. The calculation of the integral is still rather complicated. In order to compare 
the analytical solution with the numerical results, we can approximate this moving 
boundary problem as a slab with an increasing thickness with time. By employing the 


solution for heat conduction of slab with fixed thickness and replacing the thickness 


by A(t) = Ut, we have 


Swe aye 3 (I) erie Ta z jer Hh, (3.31) 


with 
O0<z< Ut. 


This is easy to calculate and the sum of the first several terms gives enough accuracy 


to compare with the numerical results. Now we consider two special cases. 


Slow conduction (A << 1) 

In this case, only the first term in the terms when n = 0 in the above solution 
is dominant. All the other term vanish very quickly. This corresponds to the semi- 
infinite space solution for the heat conduction with a constant heat flux at z = 0 


(Carslaw & Jaeger 1959). That is 
AAt z 
Q(z, t) = (——)/7ierfe ———_.. 300) 
(2,1) = (=) = (3.32) 
Fast conduction (A >> 1) 


When A >> 1, the temperature will quickly reach a steady-state. The temperature 


equation is approximately 


mle) 
—=(0 3.33 
Oz , 238) 
with 
Oo 1 
O42] 01) = and pz he = 9) ai (3.34) 
The solution for the equation can be easily obtained. We have the steady state 
solution 
Ut — 
Q(z,t) = —~*. (3.35) 


Kk 
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The comparison of the above temperature solution is given in Fig. 3.9. The dashed 
curves are the asymptotic solutions, the asterisk (*) corresponds to the data of steady 
state solution, and the solid lines are the numerical results. It is clearly seen in this 


figure that they all are consistent. 


3.4.5 Effect of diagenesis 


Effect of diagenesis on compaction: illite & montmorillonite fractions in Fig. 
3.10 show illite formation and the montmorillonite diagenesis process. Values of 
A= 1,A = 1,R = 0.01,a = 0.3,n = 10 and t = 5 have been used. The extent 
and speed of diagenesis essentially depend on the temperature and residence time of 
diagenetically active temperature. Diagenesis proceeds more efficiently and fully at 
a higher temperature and deeper burial depth than that at a lower temperature and 


shallower region. 


Effect of diagenesis 


montmorillonite 


I 


0 0.1 0.2 0.3 0.4 0.5 
Porosity 


Figure 3.10 Effect of diagenesis on compaction. Z is the scaled height. 
Dashed curve corresponds to the solution for the case of no diagenesis. 
Mechanical compaction is clearly the main important factor controlling 
the porosity evolution while diagenesis is of secondary importance with 
this choice of parameter values. 

This figure presents a more complete and full view of porosity evolution during 


diagenesis. From this figure, we see that mechanical compaction is the most important 
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factor controlling the porosity evolution, while diagenesis is also a very important 
factor, but it is in the secondary position. 

In this chapter, we have only provided some numerical results for the cases of 
geological importance to indicate some features of the compaction and diagenesis 
processes. A mathematical analysis of the model is the main purpose of the following 


chapters. 


Chapter 4 


Asymptotic Analysis and 


Comparison 


Despite the importance of compaction, few analytical solutions are available for situ- 
ations of practical importance. Gibson (1958) obtained a solution in terms of excess 
pressure with prescribed constant moving boundary. This solution is most commonly 
used in the literature for linear compaction theory. Audet & Fowler (1992) obtained 
two asymptotic solutions for the case of 4 >> 1 and t >> 1, and for the case of 
A << 1. But these solutions do not fall into the time ranges of geological interest. 
In fact, the useful solutions for the evolutionary history are those with t = O(1). 
Hence the main purpose in this chapter is to extend Audet & Fowler’s work to the 
geologically relevant situations of smaller times. 

In order to verify the validity of the numerical results, we will use an asymptotic 
analysis to give an approximate description of the different cases corresponding to 
different values of the dimensionless parameters and compare the numerical results 
with the approximate solutions with the same parameters. The analysis is mainly 
based on the relative size of A and m. The analsysis is further elaborated in Fowler 
& Yang (1997) who investigate slow and fast compaction in sedimentary basins and 


the related geological significance. 
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4.1 Non-linear Diffusion Equation for Porosity Evolution 


We now neglect diagenesis (k, = 0 or R = 0) and let a = 0, and restrict our attention 
to a single solid species. Let ¢ = ¢; be the porosity, then the porosity equation 


(3.16)-(3.19) degenerates into a general non-linear diffusion equation 


a6. 8 


fens oe L = 2 pe 
Fy 7 gs thl - HGH - UI} (4.1) 
k = (¢/do)™, m>> 1, (4.2) 

- ~ 10¢ 
ania laa) erage (4.3) 
with boundary conditions 

¢,-¢=O0at z=0, (4.4) 
o = do at z=h. (4.5) 


This is a non-linear diffusion problem with a free boundary, whose behaviour is es- 


sentially controlled by the dimensionless parameter A. 


4.2 Analysis 


From the parameter estimation, we understand that values of A will usually lie in 
the range 10-7 — 103. Since 2 is the controlling parameter which characterises the 
compaction behaviour, we can therefore expect that \ = 1 defines a transition between 
slow sedimentation 1 >> 1 and fast sedimentation 4 << 1, and that the evolution 


features of fast and slow compaction may be also quite different. 


4.2.1 Slow compaction (\ << 1) 


For A << 1 and z ~ 1, the @ equation implies that 0¢/Ot % 0, with ¢[h(t)] = do, 
therefore, ¢ © ¢y and k = 1. The outer solution ¢ & $9 does not satisfy the boundary 
condition at the base z = 0, which implies that there exists a boundary layer near 


z = 0. From the numerical results, and the fact that 4 << 1 corresponds to the fast 
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sedimentation case, we then can write the ¢ equation in the form in terms of rescaled 


inner variable ¢ 


06 Pe 2 — A (l= $0)? 


at OC?’ eG Steer do” ee) 
with boundary conditions 
oe = jue C=, (4.7) 
OC 
@—> do as Co. (4.8) 


This @ equation is equivalent to the case of heat conduction in a semi-infinite space 
with a radiation boundary at z = 0 and with a far field matching condition (equivalent 
to the initial temperature condition). The solution can be easily obtained by the 
standard Laplace transformation method (Carslaw & Jaeger 1959) 


o= overt ayia! + bret Merfel ra EMEA]. (4.9) 


This solution shows that for the case of A << 1, the sedimentation is so fast that 
the compaction can only develop in a small range near the basin basement with a 
thickness proportional to V’/t. When a = 0, we are in the case discussed by Audet 
& Fowler (1992) with a similarity solution (their equation (5.26)). 

Audet & Fowler’s solution (5.26) is in fact equivalent to the case of conduction in 
a semi-infinite space with a constant flux ¢, = ¢ at z = 0 into the medium with 
zero ‘temperature’. The solution of this case can be expressed exactly as (Carslaw & 
Jaeger 1959) 


heat. me a acs z Po 
} = bo — doV4N't ierfc(E), € = ee (4.10) 


where 


ierfe(€) = ae — €erfc(€). (4.11) 


This solution is essentially the same solution as equation (5.26) given by Audet & 
Fowler (1992). Audet & Fowler’s solution is only an approximation with a constant 
flux boundary, which is accurate if V\’t << 1. As \ << 1, we expect that this 


solution is a good approximation when t < O(1). If t is large, then this solution 
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will break down. But the solution (4.9) will hold uniformly for all time. In fact, if 
Vt << 1, Both equation (4.9) and Audet & Fowler’s equation will approximately 


predict the same value at z = 0 (i.e. € = 0) 


o(z = 0) © d0 — do4/ = (4.12) 


When € is large (€ — 00), by using the asymptotic expansion of erf(€) (Hinch 1991) 


ee 1 : 
— La +...) with €— oo, (4.13) 


22 
and € >> VA/t, we can write both solutions in the same approximate expression 
V AN t —€2 
ake 
262/74 


The comparison of the above solution with the numerical results is given in Fig. 


erf(é) = 1 


(4.14) 


~ & bo — oo 


4.1. This shows the good agreement between the solution (4.9) and the numerical 
results. The agreement between Audet & Fowler’s solution and the numerical results 
is almost the same as the solution (4.9) when t is small, but it clearly gets worse when 


t becomes larger. 


Comparison of diffusion type solution with numerical results 
T T T T T r T 
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036. 038 04 0 oad 046 048 
Figure 4.1 Comparison of analytical solutions with numerical results 
(solid) (A = 0.01). The diffusion solution (4.9) (dashed) and Audet 
& Fowler’s solution (dotted) are plotted versus the scaled height Z at 


different values of time t. 
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The comparison suggests that the mechanism of porosity evolution for the case of 
A << 1 is essentially controlled by the diffusion mechanism. In the present case, the 
problem is equivalent to the case of heat radiation into a semi-infinite space at z = 0. 


The overburden, hydrostatic and excess pore pressures satisfy, respectively, 


—altr-¢ (4.15) 
_ =i (4.16) 

ODa 
~F* = (1- 6)(1- 64/9). (4.17) 


For the case of \ << 1, we substitute the solution (4.9) into (4.17) and integrate 
from A(t) to z with a boundary condition p, = 0 at the top z = A(t), to obtain 


Pa = (1 — 0)(h — 2) — (1 — ¢)[erf(h) — erf(z)| 
h 


—(1- POOH Tana pl Nee = erle( se MEE, (4.18) 


This solution gives the leading order solution py © (1—@9)(h—z). The other terms are 


only small corrections. The excess pressure develops proportionally to basin thickness. 
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Figure 4.2 Hydrostatic, pore and overburden pressures at t = 5 for 
A = 0.01. Solid lines correspond to numerical results, the dashed line 
is calculated from solution (4.18). The numerical and analytical results 


are indistinguishable. 
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The comparison of the above solution with the numerical results is plotted in Fig. 
4.2. The pressure is rescaled in such a way that p;, = r(= 0.6) at the base. It can be 
seen that the agreement is very good, and that for \ << 1, overpressure is essentially 


proportional to basin thickness. 


4.2.2 Fast compaction (\ >> 1) 


From the above @ equation (4.1), it is clearly seen that the control parameter . is 
always combined with k. This suggests that \k = 1 will define a transition for the o) 


solutions. This condition gives equivalently a critical value of ¢* 
o* = doe ™. (4.19) 


Thus ¢ > ¢* corresponds to 4k >> 1 which is the range of z ~ h(t) at the top of 
the basin, while ¢ < ¢* corresponds to Ak << 1 which is the range near the bottom 
of the basin. The features of the solution in these two ranges can be expected to be 
different. For t less than a critical value to, there is not enough time for compaction 
to proceed, then we will have ¢ > ¢* everywhere, so that the low Ak regime will only 


exist for t > to. 


4.2.3 Compaction of thin sediment layers (¢ > ¢* with t < to) 


When 4 is large, the problem is one of singular perturbation type. We will assume 


expansions of the form 


1 1 

6 =o) ae a2 ee ae (4.20) 

p= A 4+ 2A 4 +,@ (4.21) 
i o ue 


Substituting the above expansions into (4.1), and equating the coefficients of powers 


of 1/X, we have 


O + 
5 {ha - PSG a) — 1} =0, a 
fa) = 1 (1) 
=e OY Flot? 50) ee 
es (2) 
4 = RU PTO? — THON (4.24) 
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where k® = (© /¢9)™. The boundary conditions become 


at z = AO 
go = o, 
dD 4 pMgO = 0, 
1 
go?) + see =0,..., 
ab 20 
(0) _ 40) 
od, =e, 
@y..4@) 
Qs Q 9 
2) = GO 
with 
gM 


[o? _ goer | 


1 


pO =1+k(1 — 6) 50 


on z = A), 


Integrating equation (4.22) and using boundary condition (4.26), we have 
~ 1 
ke(1— aloe —lj=0. 
Since k° 4 0, we have 


1 
Its solution is then 


p = dye PO), 
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(4.25) 


(4.26) 


(4.27) 


(4.28) 


(4.29) 


(4.30) 


(4.31) 


This leading order solution is essentially a steady state solution which corresponds 


to compaction equilibrium to which the porosity curve will tend when t — oo. This 


exponentially decreasing solution was obtained by Athy (1930) by fitting the observed 


data of Paleozoic shales from Kansas and Oklahoma. Athy’s porosity curve repre- 


sents compaction equilibrium attained over a very long time span. Hedberg’s (1936) 


porosity curve for the Tertiary shales in Venezuela is similar to Athy’s curve. 


CHAPTER 4. ASYMPTOTIC ANALYSIS AND COMPARISON 49 


From equation (4.1), we notice that the perturbation method is only valid if Ak >> 

1, ie. exp{m[II — (h© — z)]} >> 1 where II = (In \)/m. If \ = 100 and m = 8, then 

II © 0.58. Therefore, the leading term solution 6 is expected to be valid under the 
condition 

z>h© —IL (4.32) 

The comparison of the solution with related numerical results is presented in Fig. 

4.3. The comparison clearly shows that Athy’s relation (Athy, 1930) of porosity- 

burial depth is only valid in the range of 0 — 0.58d km in such sedimentary basins 


where their control parameter \ >> 1. If d = 1 km, then the range is 0 — 580 metres. 


Comparison of Athy-type Solution with Numeric Results 
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Figure 4.3 Comparison of Athy-type solutions (dashed) with numerical 
results (solid) for A = 100. Z is scaled height, and to is the time given 
by equation (4.39). 
Using the solution (4.31), equation (4.23) becomes 
Sean 
pO) 
Integrating the above equation, using boundary condition (4.26), and noticing that 


© /o© = 1, we have 


; ‘ a. 
AO dye OOP = {HC — GOP Ty AM — OM}, (4.33) 


a) hOdo(1 — e)eh® g() 
ko — 6 


oY — ¢ =) (4.34) 


Using this equation, equation (4.27) and solution (4.31), we obtain a relation for h 


= =A) =¢o) + AM Ge") =: (4.35) 
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Integrating this equation, we have 
nh = (1— do)t + doll — e™”). (4.36) 
Clearly, if ¢ is large, then exp[—h ] << 1, we thus have 
Ah = 1 - do. (4.37) 
If t is small, then exp[—h] ~ 1, we have 
AO 1. (4.38) 


Putting A = Tl, we can obtain an explicit expression for tp by using the above 
solution: 


fee me (4.39) 


If \ = 100 and m = 8, then II © 0.58, tp © 0.71. 
The solution of equation (4.34) with boundary condition (4.25) is then 


bY = doe7[—h Me” — y(0, 2) + x (0, 2)], (4.40) 
where és a 
7 (eR et een 


By using m >> 1, this integral can be approximately expressed as 


Po e-hO> m(hO—2)—pO 
(0.2) = Seiya a gy lld — mm + me Nene 
—(l—m+ rae em Re (4.42) 


Substituting this integral into (4.40), we have 


pi) = Gee + 


Po 
m(m — 1)(1 — 0)? 
By using equation (4.24), boundary condition (4.26) and solution (4.31), we can 


[(l—-—m+ mer Pm amnion — as (4.43) 


obtain a relation for h® at zx = h© 


. h (9) 
hO(1 — $9) = ob. dz. (4.44) 
0 
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Employing solution (4.43) and noticing that e~(™-)" = 0 when m >> 1, we have 
ploying g 
an approximate solution for h) 
j,(0))2 
pong jue (4.45) 
m 


From the above solution, we find that the perturbation method is only valid for 
sufficiently small t, otherwise h goes unboundedly. 

From the above solution and the numerical results, we find that for the case of 
t < to, the porosity has not reduced to a value ¢ < ¢*, so the case @ < ¢* with t < to 


need not be considered. 


4.2.4 Compaction of thick sediment layer (¢ < ¢* with t > to) 


Note that, from the definition, ¢* << 1 if A >> 1, so that we must formally assume 
m >> 1 in order to have ¢* of order one. Thus, we now consider a limit in which m 


is large. For convenience in the following discussion, we set 
pipe eles : * —1tind 
og=¢e =~ with d*=¢de™™. (4.46) 
Then the ¢ equation becomes 


0 1 
Spell” = 5, {el - gem) wb, — I}. (4.47) 


Noticing that m >> 1 and exp(q/m) = O(1), the above equation is then simplified 


as 
1 — o* 2 
with boundary conditions 
p=” at 2z=—h bt) — IL (4.49) 
and 
sy = Hie ab 0: (4.50) 


From the method of characteristics, we have 


w=0 and z=Ke?. (4.51) 
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The solution satisfying the boundary condition is then 


wy a w(T), 
z= Ke’ (t— 7), (4.52) 


where (7) will be determined by the boundary condition at z = 0. From the above 


solution and the boundary condition (4.50), we have 


Wz =U, (T)-T2 =m, (4.53) 


which is simply 
wt mK — 0. (4.54) 
Integrating this equation, and using that when 7 = 0, w = 0 and h = II which 


corresponds to the fixed time t = to, we obtain 


i 
mK(r —to) +1 


u(r) = In| iP (4.55) 


Substituting 7 from solution of (4.52) into the above solution and rearranging the 


equation, we have 
1+mz 


trl ; 4.56 
ved) = ble (4.56) 
Using (4.46), we finally have 
Ae Lh haa 
Zi SO: - -—|m, ALOT 
Wet) =r gal ie (4.57) 


The fixed time ty, which is given by equation (4.39), defines a lower time value under 
which the solution will be invalid. 
When ¢ is large (i.e. t >> tg, z = O(t) >> 1), then the solution (4.57) can be 


expressed approximately as 


b= o'()m, (4.58) 


pz 
mt 


Using the definition of ¢* in (4.19) in the above expression, (1/\)/™ ~ (1/A)/"—D 


as m >> 1, and putting € = z/t, we have 


b = oo(—) =. (4.59) 


This is exactly the same solution obtained by Audet & Fowler (1992, equation (5.9)) 


for the case of A >> 1 andt >> 1. 
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4.2.5 Compaction of thick sediment layers (¢ > ¢* with t > to) 


In this case, the equation (4.29) will not be valid, and a more general expression is 
1 
k(1 — WV Foe —l] = F(t), (4.60) 


where F(t) is a function of t only. From the moving boundary condition, we have 
F(t) = (h—-1)(1— ¢p). (4.61) 


This is only valid in the region with a depth less than II from the top boundary. In 
the region which includes the transition region of @ * ¢*, the term ¢; can not be 
ignored. The three terms in the ¢ equation must be considered at the same time. In 
fact, from the leading solution (4.31) in the perturbation method, we have ¢; ~ —hd,. 

From (4.31), we find that ¢@ depends on h(t) — z near the top, ie. 6 ~ (h(t) — z). 
From the numerical results, we observed that @ decreases nearly exponentially with 
increasing depth 7 = h(t)—z in the top region. This suggests a solution for ¢ equation 
in the form 


¢=¢(n), with n= h(t) —z, (4.62) 


then the ¢ equation (4.1) becomes 


I wt m 1 ! f 
hdl = (Lyn a - eGo +0 (4.63) 
Po ? 
where a prime means a differentiation with respect to 7. The boundary conditions 
(on 7 = 0) are 
@ = 0; 
7 ) m 1 ! 
h=1-X(—)"(01 - 6)(-—¢ +1). (4.64) 
Po ) 


We can see that (2.64) will imply h =const due to (4.31) and (4.62). By integrating 


the above equation again and using its top boundary condition at 7 = 0, we have 


hbo — #) = (1—A)(1 — bo) — “yn = Ge +1), (4.65) 


whose solution can be written as a quadrature. The undetermined h in this solution 


will be determined by matching it to that in a transition layer analysed below. 
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The comparisons of the travelling wave solution (4.65) and solution (4.57) (dashed) 


with the numerical results (solid) are shown in Fig. 4.4 and Fig.4.5 


Comparison of Asymptotic Solution with Numeric Results 
T T T 


1 T 
R=0, a=0, lambda=100 
0.9 
Dashed: Solution 
0.8 
Solid: Numeric 
0.7 
it 
1 t 
0.6 ay a 
! ' - 
1 
N 0.5 i 
if 
ii 
0.4 iff 
if 
ift 
0.3 ft 
i 
f 
0.27 J ; 
Pa; 
0.1 /| 
y 
t=5,/ 7/ 2 
0) 1 raf ad 1 1 1 
0 0.1 0.2 0.3 0.4 0.5 


Porosity 


Figure 4.4 Comparison of asymptotic solutions (dashed) with numerical 
results (solid) at t = 2, 5 for \ = 100. The upper two dashed curves 
correspond to solution (4.65), and the lower two correspond to solution 


(4.57). Z is scaled height. Agreement gets better as t increases. 
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Figure 4.5 Comparison of asymptotic solutions (dashed) with numerical 
results (solid) with different values of m. The dashed curves have the 


same meaning as in Fig. 4.4. The agreement gets better as m increases. 


The porosity profile for t = 2,5 is plotted in figure 4.4 with other values fixed 
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(a = 0, A = 100, R = 0). In Fig. 4.5, different values of m are compared for the 
same t = 2 with all the other values fixed as in Fig. 4.4. They clearly show that their 
consistency becomes better for larger times or larger m even with small t. This is 


what we have expected from the condition used to derive the solutions. 


4.2.6 Matching the solutions 


In order to match the solution (4.65) to the solution (4.57), we define a transition 


region by adopting the transition variable ¢, 
1 
zg=A(t) -II+ = ie. ¢ = m(h(t) — II — 2). (4.66) 


Rewriting the solution (4.57) in terms of the new variable with UV = w as in (4.46), 
we have 


g=¢'e mm, (4.67) 


where 
1+m(h—II)+¢ 


UV =I1n . : 4.68 
ae Ftp) +1 ae 
Noticing that m >> 1, we have approximately, for the lower solution (4.57), 
h —II 
teh) (4.69) 
G 2 (t = to) 


for 1 <<< —¢ << m. Now the ¢ equation in the transition region can be written as 


—mhgt [et] = moe (1 — $*)*(Ue ~ 1) (4.70) 


Integrating this equation, we find 


v-lInm 


~e¥(1 = 6° )*(We = 1) + hgte 


Woo—Inm 
m 


—hd*e - + e¥~(1 — g*)?, (4.71) 
m 


where we require V — UV, as ¢ — —oo. Comparing with (4.69), we have 


—Il 
ge (bt) 


or 


Voo-—Inm 


ee ae ee (4.73) 
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Rewriting and rearranging the upper solution (4.65) in terms of the new variables 
¢ and W, we have 


w-—lInm 


~¥(1 — $)(Ue 1) + hdte B® = —(1 — bo) $A (4.74) 


By using Van Dyke’s matching rule (Van Dyke, 1964; Hinch, 1991), we expect that 
the W obtained from (4.71) and (4.74) should be the same in the matching region. 
From equations (4.71) and (4.74), we notice that the left hand sides of both equations 
are the same and independent of ¢, thus we have 


co—Inm 


—(1—¢))+h= hdte~ mm * = (1 -—¢y. (4.75) 


using the fact that m >> 1, and rearranging (4.75), we finally obtain an equation for 
h(t): 


i 

pe odo 
i Poe 

which determines A(t). It is worth pointing out that (Y —Inm)/m in the second term 


(4.76) 


of the left hand side of (4.74) is not accurately set to zero, since it is order of — 4 Inm. 
But if we do set (V —Inm)/m ®& 0, then we obtain the leading order approximation 


for h: 

a VLOG 

he : 
1-¢ 


Clearly, the non-negligible term (WY — Inm)/m will provide us a more accurate ap- 


(4.77) 


proximation for h. 

Now we understand that the top solution breaks down as 7 > II while the bottom 
solution fails as z > h(t) — I]. We can simply construct a uniformly valid asymptotic 
solution (Hinch 1991) since the solution in the upper region is the same as that in 
the lower transition region. If we note the solution in the top region as ¢.) and that 


in the bottom region @pottom, the composite solution is then 


o(z, t) = Dbottoni + Ptop aa Deo: (4.78) 


From equations (4.59) and (4.73), we know that @,, is time independent when t is 


large (t — 00), therefore, h is a constant. We simply have 


ie : =O with Ox, © dul Oe. (4.79) 
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Clearly, this is exactly the equation (5.16) obtained by Audet & Fowler (1992). In 
this case, the solution is a travelling wave solution. 
The comparison of the matched composite asymptotic solution (4.78) with numer- 


ical results (dashed lines) is shown in Fig. 4.6. 
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Figure 4.6 Comparison of matched asymptotic solutions (dashed) with 


numerical results (solid) at t = 2, 5 for A = 100. 
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Figure 4.7 Comparison of h(t) solutions (dashed) with numerical results 
(solid) with \ = 100. 


For the case of X >> 1 and ¢ > ¢*, substituting Athy’s solution for ¢ in the top 
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region, we have 


Pa 

> i; (4.80) 
This equation with the boundary condition p, = 0 at the top z = h(t) gives that 
Da = 0 at the leading order. This means excess pressure does not occur for short 
times or in the top region where z > h(t) — II. This region is clearly shown in Fig. 


4.8. For larger times, the solution suggests that 6, << $, whence 


ODa = 
Pe ws (1-4), (4.81) 


which shows that the excess pore pressure develops at large times even if \ >> 1. 
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Figure 4.8 Hydrostatic, pore and overburden pressures at t = 5 for 
A = 100. Dashed curves are computed by using (4.65) and (4.57). 
The comparison of the numerical results with the pore pressure calculated from the 
asymptotic solutions (dashed) is shown in Fig. 4.8. The overpressure only develops in 
the lower region, while the pore pressure remains hydrostatic in the top region with 


a depth of order I] from the surface. 


4.3 Summary 


In summary, we find that the limit \ << 1 (slow compaction) can be simply analysed 
by means of a boundary layer analysis at the sediment base. Essentially, sediment is 


added so fast that the porosity remains virgin except near the base, where compaction 
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occurs. The pore pressure is then essentially lithostatic, that is, excess pore pressures 
exist over the whole domain. 

The more interesting mathematical case is when A >> 1 (fast compaction). For 
sufficiently small times, the porosity profile is exponential with depth, corresponding 
to an equilibrium (long-time) profile. However, because of the large exponent m in 
the permeability law k = (6/0), we find that even if \ >> 1, the product \k may 
become small at sufficiently large depths. In this case, the porosity profile consists of 
an upper part near the surface where \k >> 1 and the equilibrium is attained, and a 
lower part where \k << 1, and the porosity is higher than equilibrium. Straightfor- 
ward asymptotic methods are difficult to implement because the limit m >> 1 implies 
exponential asymptotics, but we use a hybrid method which appears to correspond 


accurately to numerical computations. 
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Fig. 4.9 Regions of different solutions (A >> 1). A, D, E are the 
regions with @ > ¢* while B, C are the ones with ¢@ < ¢*. The region 
F between the two dashed lines is the transition region with @ = @*. 
Audet & Fowler’s regions C, E are for large times (t — 00). 
To summarise the solutions for the case of \ >> 1, we represent the solutions in 
their related valid regions in Fig. 4.9. The regions below h(t) line labelled as A, D, E 
are the regions with ¢@ > ¢* while those labelled as B, C are the ones with ¢@ < ¢*. The 


region with dashed lines on both sides is the one with ¢ & ¢* which is the transition 
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region. The regions E, C on the right side of the dotted line constitutes Audet & 
Fowler’s region (1992) which is only valid for large times (t — oo). 

Correspondingly, the terms in the ¢ equation will play different roles in the be- 
haviour of the solutions. The over- and under- braces label the relative terms which 
govern the features of the solutions in different regions. 


RegionA,D&E 
ee N 


06.8. » , , 106 
ae er pees 82 
Fe = ge tRl— 9)" 14+ 552) (4.82) 
RegionB&C 
When the left hand side ¢; is negligible, we have 
i 2 106 
- = a Peer U. 4, 
{RL 9)"[-14 SEI}. ~0 (4.83) 


then we have the solutions for the top regions (¢ > ¢*). If the first integral of the 
right hand side is zero, the Athy-type solution is obtained in region A. If the first 
integral is not negligible, the solution in region D is thus obtained. If ¢ is large, this 
solution moves into the region E of the travelling wave type solution which is given 
by Audet & Fowler (1992). 

When the diffusion term on the right hand side is negligible, we have 


ob, & —Mk(1— ¢)?}., (4.84) 


thus the solution for the bottom region (@ < ¢*) is obtained. The limit for large t 
of this solution is exactly Audet & Fowler’s solution for large times. In the region 
F (¢ & @*), all the three terms in the @ equation must be considered. The matched 
composite asymptotic solutions provide a uniformly balanced solution for the whole 
region. 

The methods presented in this paper pave the path for the analysis of compaction 
in sedimentary basins when more complicated loading histories are studied, and also 
when more realistic phenomena are included, such as diagenesis, or state-dependent 


rheology. 


Chapter 5 


Unloading and Variation of 


Sedimentation Rate 


In the model we analysed in the previous chapters, the rheology of the porous medium 
is considered as poro-elastic, and it is equivalent to a single-valued function of the 
Athy’s type pe = pe(@) in the 1-D case. A more realistic rheological relation should 
include the nonlinear effect of hysteresis derived from soil tests. In addition, the 
sedimentation rate m, has also been taken as a constant in the poroelastic compaction 
model. From the numerical simulations in chapter 3 and the analysis in chapter 4, 
we can see that the model does not require the sedimentation rate to be constant. In 
fact, the dimensionless sedimentation rate m, can vary with time t and it can also be 
negative, which corresponds to the case of unloading. 

In this chapter, we will mainly investigate the effect of unloading and variation of 
the sedimentation rate by using a more realistic rheological relation. As the analysis 
of the model equations is very complicated, we will simply show the numerical results 


and give some analysis whenever possible. 
5.1 Model Equations for Unloading and Reloading 


5.1.1 Non-linear soil behaviour 


In order to model the phenomena of unloading and reloading, we must consider 


the non-linear stress-strain behaviour which has been investigated in many cases. 
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Bethke & Corbet (1988) examined the non-linear effects associated with the long- 
term compaction of sedimentary basins. The one-dimensional isotropic consolidation 
test (Das 1983, Burland 1990) of soils accompanied by unloading/reloading sequences 
clearly shows that the soil behaviour is path-dependent and nonlinear as shown in 
Fig. 5.1. The behaviour during unloading and reloading is essentially elastic with a 


small amount of hysteresis. The void ratio 


e=4/(1-4) (5.1) 


is used in this figure as the conventional way of presenting the test results. 

In order to model the behaviour of soils as shown in this figure, the Cam-clay mod- 
els developed by the Cambridge group, in terms of Critical State Formulations, are 
very attractive since these models are able to reproduce qualitatively a good number of 
the main features of the mechanical behaviour of soils such as unloading/reloading, 
stress path-dependence etc (Schofield & Wroth, 1968; Atkinson & Bransby, 1978, 
Huekel & Baldi 1990). If a more accurate reproduction of actual soil behaviour is 
sought, the more sophisticated models such as the Modified Cam-clay model (Roscoe 
& Burland, 1968) and the more modern cap model (Chen & Mizuno, 1990) should 


be used. 


v=l+e v=l+e 


Inn Ing 


(a) (b) 
Figure 5.1 Non-linear behaviour of soil consolidation. NCL is the nor- 


mal consolidation line, and URL is the unloading and reloading line. 
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The modified Cam-clay model is an isotropic, nonlinear elastic strain-hardening 
plastic model in which only volumetric strain is assumed to be partially recoverable. 
The consolidation curve in Fig. 5.1(a) is idealised as that shown in Fig. 5.1(b). The 
virgin isotropic normal consolidation line is assumed to be linear, and the unloading 
and reloading curves are parallel as a single straight line. 


The equation for the normal consolidation line (NCL line) is given as 


€ = €) — C_In(p-/Do), (5.2) 


where C, is the compression index (Das, 1983; Burland, 1990). For the unloading- 
reloading line (URL line), we have 


e = €) — C,ln(pe/po), (5.3) 


which is valid when pe < p*, where p? is the maximum of previous values of pe. 
Similarly, C, is called the swelling index. The join of the two lines corresponds to 
a special value p* of p. in the time-history of compaction, and the value ej can be 


expressed in terms of p= as 
€ = €0 — (Ce — Cs)In(p2/po). (5.4) 
If written in dimensionless form in terms of p, then equations (5.2) and (5.3) become 


€ = €) — C_In(p/fo), (5.5) 

and 
e = 9 — (Ce. — C5)In(p" /Po) — Csln(p/Po). (5.6) 
In equations (5.5) and (5.6), we provide only one of the many possible formulations 
of the nonliear constitutive laws which can be derived from the modified Cam-clay 
model (Roscoe & Burland, 1968) and modern cap model (Chen & Mizuno, 1990). In 
fact, what we have used before in equation (3.1) is just another form of the formula- 
tions, and is widely used in the literature (Smith, 1971; Sharp, 1976; Das, 1983; Audet 
& Fowler, 1992; Wangen, 1992). Different formulations will result different forms of 
function p(@) or e(p), but they all can reproduce the main features of the nonlin- 
ear behavior of loading/unloading. In order to compare with the results obtained in 


Chapter 4, we will use the modified constitutive law (5.15) similar to (3.1). 
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5.1.2 1-D model equations 


From the derivation of the general model discussed before, we know that the model 
does not require the sedimentation rate to be constant. Nor does it require increasing 
loading only. Change of sedimentation rate and erosional unloading can be treated 
within this model, but the constitutive laws for the case of unloading should be 
changed. 

Recalling the process of non-dimensionalization in Chapter 2, we defined a length- 
scale d in (2.48), and scaled z with d, u*® with m,, time t with d/m,, pore pressure p 
with (9,;—pz)gd, and permeability & with ko. Clearly, the scalings involve the sedimen- 
tation rate m,, and thus must be modified to allow the variation of the sedimentation 
rate and erosional unloading. 

If we scale time t with a time scale 7 instead of d/m,, u® and m, with d/r instead 
of m,, and keep all the scalings of other quantities the same as before on Page 22, 
then the obtained dimensionless governing equations are the same as equations (2.51)- 
(2.58). The only change is to replace the m, by d/r in the expressions of A, A and 
R. Thus (2.59)-(2.61) are replaced by 


y — holes = pig Ko knd 


udjr) ? peddfp@’ ~~ (afr) 


Clearly, if one substitutes mm, = d/7 back into the above expressions of A, A and 


(5.7) 


R, we do have the same expressions as (2.59)-(2.61) in the case of constant sedi- 
mentation rate. Therefore, the dimensionless model equations (2.51)-(2.58) are still 
suitable when sedimentation rate changes, but the real meaning of m, is the aver- 
age sedimentation rate in the relevant time history of sedimentation. In the case of 
constant sedimentation rate, m, is the real constant sedimentation rate. 

In order to show more efficiently the effect of the variation of sedimentation rates 


and erosional unloading on the porosity evolution, it is convenient to ignore diagenesis 


and temperature effects by setting R = 0, a = 0, b = 0, ¢ = ¢ in equations (2.54) 


and (2.58), and omitting the temperature equation (2.57). By using the force balance 
equation (2.58) to eliminate p in Darcy’s law (2.56), and using the expressions (2.64)- 
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(2.65) for u’ and u!, we finally obtain a single non-linear diffusion equation 


06 .0.; Op 
se = AS {ha - HI-F - 1-9), (6.8) 
with constitutive laws 
p=p(d), k=k(¢), (5.9) 


which are given below in equation (5.13)-(5.15). 


Boundary conditions 


Rewriting the definition of the effective pressure in dimensionless form, we have 
p=P-p. (5.10) 


In the case of very rapid unloading (Haxby & Turcotte, 1976), P decreases suddenly, 
but p may not have enough time to respond to such a quick change, and thus remains 
nearly a constant, which subsequently forces the effective pressure p < 0. The whole 
column of the sediments will be unloaded instantaneously. The negative effective 
pressure implies that fracturing should occur, and the model equations will become 
invalid for fracturing. In reality, the unloading due to erosion at basin surface is a 
very slow process, and the effective pressure should be always non-negative, p > 0. 
Therefore, a reasonable boundary condition at the basin top z = A(t) in the present 
model is to assume that the effective pressure p always remains zero, i.e., p = 0, 
which eliminates the possibility of fracturing due to very quick unloading, discussed 
by Haxby & Turcotte (1976). 


Now the boundary conditions are 


—5, ~U- 4) =Oat z=0, (5.11) 
og=¢o at z=h, 
h = m(t) + yee —(1-4)| (5.12) 
Oz , 


Here, 7m is the dimensionless sedimentation rate which is 1 if it is constant, or O(1) 
if time-varying. It is based on these equations that the change of sedimentation rate 
(m(t) > 0, increasing loading) and erosional unloading (m(t) < 0) will be treated, 


but the constitutive laws will change correspondingly. 
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5.1.3 A specific case 


To investigate the main features for the cases of our interest and compare with the 
earlier results in Chapter 4, we still use the Smith type constitutive function of per- 
meability as before, 


k = ($0/¢)" with m=8. (5.13) 


But the constitutive relation (3.1) 


P = In(¢0/¢) — (¢0 — 4), (5.14) 


is only valid on the increasing loading branch (NCL). On the unloading-reloading 


branch (URL), we use the following constitutive relation: 


POY 


pale V) Gola ay ga 


]—(Go-@) +B", (5.15) 


with 
go(z,t) = mind(z,7 <t), p*(z,t) = maxp(z,7 < ¢), 


Cs 
B* = In(0/¢3) — (0 — 8) and y= S, (5.16) 


where 7¥ is the slope ratio of the URL line to the NCL line. The normal ratio for soils 
is y © 0.1 ~ 0.25 (Das, 1983). Clearly, equation (5.15) degenerates into equation 
(5.14) when y = 1 which corresponds to the case that URL branch falls onto, as we 
expected, the NCL branch. In this case, the behaviour of unloading and reloading is 
reversible. 

A relation similar to equation (5.15) was used by Wangen (1994), which can be 
written as 


@ = Pmin|1 aF Ole Detran — De)|; (5.17) 


but Wangen’s relation is only valid in the case of p = p* and @ & @. In fact, Wangen’s 
relation is only a special case of our relation (5.15) when y << 1. Equation (5.15) 


can be written as 
p* —p+(d—$6) 
os (1—7)¢6 ste 


ai (5.18) 
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Taking y to be small, this implies that @ & $5, whence 


yk 


Dp 
b— % 195 ae a er : (5.19) 
5 
Rearranging this equation, we finally have 
¢ = doll + — (Ph — B), (5.20) 


“% 
which is similar to the equation (10) used by Wangen (1994). 

The switch conditions for loading and unloading at any point following the material 
are 


On URL branch: 


dp* aes dp 
— f = - —_ 
a 0, if p=p and it, <Q), 
ES, Hea a (5.21) 
dt. areata Pp P ’ 
On NCL branch: 
| ae) a a ae 8 
-— if g-p —= 5.22 
ie a ee a 22) 


where d/dt, = O0/Ot + u°0/0z, and p*(z*,t) = maxp(z,T < t), where z* is a La 
grangian spatial coordinate which is related to z by dz/dt = u*® with z = z* at t = 0. 
If we use |u* |<< 1 as an approximation, then z* & z, and the material derivatives 
n (5.21) and (5.22) can be taken to be partial time derivatives. Now we have 
On URL branch: 
p=p if p<0, 
or p< Dp’; (5.23) 


On NCL branch: 
p=p and p>0. (5.24) 


If the constitutive relations (5.14) and (5.15) are plotted in semilogarithmic coor- 
dinates, we have the curves in Fig. 5.2 which are similar to the NCL and URL lines in 
Fig. 5.1. This means that the Athy-type relations in the present model are suitable 


and reasonable in reproducing the main features of soil behaviour. 
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Figure 5.2 Constitutive laws plotted in semilogarithmic coordinates. 


The curves are similar to the NCL and URL lines as in Figure 5.1. 
5.2 Numerical Method 


The numerical method in Chapter 3 is only robust to solve the model equations when 
@ and its first derivatives ¢, and ¢, are continuous. But for the present case, the 
non-linear history-dependent property of the porosity function may imply that @ or 
its first derivatives are discontinuous at the interface between swelling region (where 
@ or e increases) and compressing region (where @ or e decreases). Therefore, we 
should first ensure that the numerical method can work well in these cases. Special 


modification at the interface is needed. 


5.2.1 Finite difference implementation 


For convenience in the discussion of numerical method, the general 1-D model equa- 
tions can be simplified without losing its main features by leaving out the second 
term on the right side of the equation (5.8), so that we have 

0d, 90 Ob 


20 = 2 (oy, (5.25) 
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D(¢) = —K(1 — 9)6'(9) (5.26) 
Generally speaking, D(@) will take different forms on the NCL [| D*(¢)] and URL [ 
D~(@)] branches, and may be discontinuous at their interface. 

The erosional unloading or change of sedimentation rate at the top of the basin will 
usually generate a series of interfaces, which separate the swelling and compressing 
regions, travelling at different velocities down to the bottom. The advancing interface 
is determined by solving a compatibility equation which is usually derived from an 
integral formulation of the non-linear diffusion equation while the smooth solution 
away from the interface is treated with a standard finite-difference method. Therefore, 
the integral form of the conservation law gives a contour integral formulation along 


the moving boundary (interface) I, 


[, 1) A[D()¢.Jat + [o]az} =o, (5.27) 
where f(z,t) is any continuously differentiable function of z and t that vanishes on 
the boundary of the solution domain. The notation [ |] means [¢.] = @f — @7. Since 
the above relation is true for any arbitrary f(z,t), the integrand must vanish, and we 
thus obtain 

a(t) = AID()b=)] (5.28) 


This condition defines the travelling speed $(t 


_ALD( 
| 


of the interface in terms of the values 
of the solution on either side. 

To illustrate the modification of the finite difference formulae near the moving 
interface, we consider the case with a swelling region above the interface and com- 
pressing region below the interface. By using a fixed finite-difference grid, the moving 
interface, at any time j dt, will usually be located between two neighbouring grid 
points, say 7dz and (i+1)6. The unequal space intervals are used to modify the 
related finite-difference formulae near the moving interface. By using the three-point 
interpolation formulae of Lagrangian type (Crank 1975) with three known values 


f (zo), f(41), f(z2) at three points z = zo, 21, 22 respectively, we have 


fe) = oD Ly(z)f (Ze), (5.29) 
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with 
L,(z) = se ee ees 
\ 7 (z aa ZR) M4(Z)’ i 7 k=0 = 
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- ue ae 2 ee OY a ee ns = st) > wr <a} 
1 
i-l 
0z 
1-2 


Figure 5.3 Grid lines near the moving interface z = s(t). 


The formulae for space derivatives are then 


df - 2 ! fo ; ale, 
ce — 2 Ent); Li, > oy (z _— Ze )(z eS 2)T(z)’ 
and @f 2 ms 1 
ay _ 
i 2 2F =) LT nee Pe 


70 


(5.30) 


(5.31) 


(5.32) 


Applying the above formulae for the grid lines (¢ — 1) dz, 1dz and the moving 


interface z = s(t) (Fig. 5.3), we have (for z < s(t)) 


Oo = 2 Pi-1 oF Ps 


Oz? ee ag a , (i=g)Q-9” z= tz, 
and 
06 1,(l-qg)di-n @2-aQ¢i , @-24)¢s eae 
Be oe 2-q l-q 'G-peQ-9” (é) — 0. 


For z > s(t) we have similarly 


2 bs _ din Outs 
O22 (dz)?*g(q+1) qs qt 


), z=(i4+1)éz, 


and 
06 1 a+ Vos | (a+ Dos 
Oz 62 q(q + 1) q 
_1Pis2) z= s(t)+0. 


q+1 


(5.33) 


(5.34) 


(5.35) 


(5.36) 
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The above modified formulae for the space derivatives near the moving interface 
are used for the points idz, s(t) and (i + 1)éz. These formulae together with the 
usual equal space interval formulae for other points can be applied for the whole 
region at any time. In addition, the switch conditions (5.23) and (5.24) are checked 
at every point near the interface at the beginning of each time step to make sure 
that the correct branch of the constitutive relation p(@) is used for the numerical 


implementation. 


5.2.2 A test case 


To test the above finite difference formulae, it is convenient to investigate first a 
simplified non-linear diffusion with discontinuous diffusion coefficient. To compare 
the numerical results with some available analytical solution, a very special semi- 
infinite case (Crank 1975) is solved numerically. Written in the variable 7 increasing 


downward with the origin at the top, the equation is 


2 

06 _ 0% 

Ot On? 

where D = D, =constant if ¢ > ¢,=constant, D = D2 =constant if @ < ¢@,. The 


(5.37) 


boundary conditions are 
$(n = 0) = go and $(7 = 00) = boo. (5.38) 
As we mentioned before, the condition at the interface 7 = s(t) is 
s(t)— _ s(t)+ 
DiGi” = Dod. (5.39) 


Crank (1975) obtained an analytical solution for this problem 


n 


= Aerf 
P= oot a5 Dt 


0<7 < s(t), (5.40) 


and 
1) 
> s(t 


b = bo + Berfe s(t) =avit, (5.41) 


where q@ is determined by 


(ds — bo) VD, : (bs — b00)V D2 


a2/4D, a a2 /4D2 a 
€ erf wD, Oe erie, TDs 


= (5.42) 
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Figure 5.4 The advancing interface with time and comparison of Crank 
solution (solid) with numerical solutions (dashed) at t = 1. 
Taking the values of D, = 5, Dp = 1, ¢9 = 0.5, @, = 0.4, and dy = 0.1, solving 
this problem numerically, we have the advancing interface in Fig. 5.4. 
The comparison of the numerical solution with Crank’s solution is shown in Fig. 


5.4 which clearly shows that the numerical method is robust. 


5.3 Irreversible unloading and reloading 


By using the modified numerical formulae, we can study the case of irreversible un- 
loading and reloading. When calR = 0 and b = 0, equation (2.65) is equivalent to 


u’ = —d(u! — u’). Now we can write the conditions at the interface as follows 


[6] =0 and [o(u! —u*)|}=0 (ie. [u*] = 0), (5.44) 


where [f] = 0 is the physical condition of continuous effective pressure, and [¢(u! — 


u®)| = 0 is derived from the condition of no fluid stored (mass conservation) at the 
interface. The condition |p| = 0 does not necessarily imply that [¢] = 0. From Fig. 
5.1(b), we understand that [p] = 0 is equivalent to [¢] = 0 only at the interface of 
loading and unloading, but it is generally equivalent to [¢]| # 0 at the interface of 
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reloading and loading although the jump [¢] may disappear when loading proceeds to 
the extent where the two branches join again. This will be also illustrated later in 
Fig. 5.9. 

In the case of [4] = 0, [u’] = 0 or [Ak(—2 — 1+ ¢] = 0 implies pC’ = pUFe, 
From Fig. 5.1(b) or Fig. 5.2, we have p3@"/py"" = 7 or 678" /¢3™ = 7. Then, we 


have a jump condition at the interface 


ibe] = (1 — oN! = (- — 1) gURE, (5.45) 


oN“ is the value of ¢, along the NCL line while being compacted. Thus ¢, will 


where 
be continuous if y = 1. The discontinuity is a property of the irreversible compaction, 


which will be illustrated later in Fig. 5.7 and Fig. 5.8. 


5.3.1 Slow compaction \ << 1 


From the numerical results and discussion in Chapter 4, we understand that the 
behaviour of small case is relatively simple. To study its main features, we use a 
step function of constant loading and constant unloading. The numerical results are 


shown in Fig. 5.5 at different times. 


Loading to t=5, unloading to t=6 


0.67 


0.57 


height (z) 


0.4 0.42 0.44 A F 0.46 0.48 0.5 

‘orosity 
Figure 5.5 Porosity profile under constant unloading (A = 0.01, 7 = 
0:25).; 7 Sa bs oa eo a Oe SL a 6 Ty 


Unloading begins at t = 5. z is the height measured from the basement. 
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We see that the unloading at the top has very little influence on the porosity profile 
in the whole region. This is actually consistent with the boundary layer phenomenon 
near the bottom we obtained before in Chapter 4. The condition for the diffusion 
boundary layer remains unchanged in the case of constant unloading and cyclic load- 
ing. The porosity in the top region is ¢ = ¢9, which lies outside the boundary layer, 
that the effective pressure there is zero; thus the pore pressure p is equal to the over- 
burden pressure P, i.e., p = P, outside the boundary layer. The change in P will 
affect the change of p instantaneously, and their changes are in phase. The effective 
pressure is only positive within a thin boundary layer near the base where compaction 
proceeds very slowly. Thus the change of sedimentation rate or unloading without 
changing its surface porosity ¢@o will not change the behaviour of the boundary layer 


near the bottom. 


5.3.2 Fast compaction \ >> 1 


Based on the previous discussion, we understand that it is the case of \ >> 1 that 
is more complicated and of more interest and importance. In order to show the main 


features of unloading, the following simple cases are investigated. 


5.3.3 Constant loading, evolving to equilibrium, then constant unloading 


We first investigate the system behaviour subject to unloading from the state of 
equilibrium. We load the system with a constant sedimentation rate to time t = 5, 
then let it evolve to its equilibrium. We then unload the system from this equilibrium 
state, and shift the time origin to t = 0 when unloading begins. The numerical results 
are shown in Fig. 5.6 in which the dashed curve corresponds to the equilibrium state. 

It is clearly seen that the interface travels downward, and is smoothed by the 
diffusion effect which only becomes important on a length scale of O(+) or O(0.1). 
Although the constitutive relations of effective pressure on porosity are two differ- 
ent functions in the unloading branch (above the interface) and compression branch 


(below the interface and at equilibrium in the present case), the model equation is 


still a single nonlinear diffusion equation whose diffusion coefficient strongly depends 
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on the porosity itself, which equivalently makes the model equation degenerate in 
such a way that the diffusive front essentially spreads with a finite speed. Thus, the 


influence region is mainly located in the top part of the basin. 


Porosity profile Interface 
3 1 
8 
2.5} £ 0.8 
So So6 
rs 5 
B15 = 
2 8.0.4 
1; , 8 
0.5: i Solid: unloading g 0.2 * data: numeric 
Dashed: equlibrium S Solid: solution 
ou 0 
0 0.2 0.4 0 0.5 1 1.5 2 
Porosity time 


Figure 5.6 Travelling interface of unloading for the case of A = 100, 7 = 
0.25. The left figure shows the travelling interface due to unloading 
at different times after unloading begins. The right figure gives the 
comparison of the solution (5.56) (solid) with numerical results (points 
with x). 
To understand this phenomenon, let us make a small perturbation w, which is 
valid at least when t is small, from the equilibrium state ¢°. The equilibrium solution 


for the @ equation when \ >> 1 is essentially the Athy type solution 
b = doe. (5.46) 


Written in terms of the (depth) variable 7 = h — z with its origin at the top, this 
equilibrium solution is 

oe = doe ". (5.47) 
Setting @ = $°+w and using equations for equilibrium state ¢* = 0, 1/|o-—(1-7) ¢°] © 
1/(y@°) since @ & ¢°, the linearised perturbation equation for w is 


OY _, 9; - 4°) oo 
a re [k re By (5.48) 
Using (5.47), we find, 
OY, 2pm) yay y= td) 
a a [e a Bn! with A=A age (5.49) 
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Le. 
ay = NeW aban — (rm — 1) ry], (5.50) 
which is only valid when t < O($). Since m >> 1, we thus have approximately 


dy + A(m — Ie" "YY, = 0. (5.51) 


For a semi-infinite space approximation, we have the initial and boundary conditions 
for w 
vin = 0,t) = f(t) and Y(t = 0,7) = 0. (5.52) 


By using the method of characteristics, we have 
w=0 and 7 =A(m—1)e""", (5.53) 


Integrating the above equations and using the initial and boundary conditions, we 


have 
it 
v= f(r) and —— F [eo _ 1] = A(m— 1)(t— 7). (5.54) 
Eliminating 7, we obtain the solution of (5.51) 
lie 5.55 
t) = f\t — ———_.5]. : 
won = fle Gaal (5.55) 


Since the unloading begins from equilibrium state, the interface travelling downward 
is the interface where w = 0 or t — (e("-)” — 1)/A(m — 1)? = 0. From the above 
solution, we therefore find that the interface s(t) is given by 


a -In{A(m ~ 1)?#+ 1} (5.56) 


m— 


Written in terms of A, we have 


t+). (5.57) 


This solution implies that the interface will travel faster as y gets smaller. In the 
extreme case when y = 0, the interface travels downward nearly instantaneously to 
the base of the column, and the whole region is unloaded (s(t) — oo as y — 0). It 
is worth pointing out that s(t) is actually the characteristics of equation (5.51) and 


thus the interface velocity s(t) is independent of unloading rate. Furthermore S(t) is 
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decreasing with time t, which implies that the unloading effect is essentially located 
in the top region at least for a short time after unloading begins. 
The comparison of this solution with numerical results is shown in the right in Fig. 


5.6. The consistency verifies the above obtained solution. 


5.3.4 Constant loading, then constant unloading 


Figure 5.7 shows the porosity profile for constant loading to t = 5, then constant 
unloading for some short times. A = 100 and y = 0.25 are fixed throughout the 
computations. 

Figure 5.7 clearly shows that an interface of discontinuous ¢, will be generated 
at the time when increasing loading switches to erosional unloading. The travelling 
velocity of the interface is not a constant. The downward travelling interface of the 
unloading region will extend the unloading region much deeper, and finally to the 
whole domain. From the numerical results, we have NC” x 0.126, 2?" = 0.032, 


[bz] © 0.094 = (1—y) oN“, which confirms that the jump condition (5.45) is satisfied. 
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Figure 5.7 The advancing interfaces of discontinuity of @, with a Heav- 
iside step function of sedimentation/erosion rate (m(t) = 1 if t < 5, 
m(t) = —lift > 5). The values of A = 100, y = 0.25 are used. Dashed 
parts shows the swelling region while the solid ones correspond to the 
compressing region at different times t = 5.1, 5.5 (or 0.1, 0.5 after un- 
loading). The right figure is the enlarged part of the left one near the 


first interface. 
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5.3.5 Cyclic loading and unloading 


To investigate the main features of the system under cyclic loading and unloading, 
a square wave function of sedimentation/erosion rates is used. Firstly, the system 
is constantly loaded to t = 5, then it goes under a square wave of unloading and 
reloading with a period T' = 1. An interesting feature arises, a discontinuous porosity 
profile as shown in Fig. 5.8 (only the first cycle is shown). There exist two travelling 


interfaces after a cycle of loading-unloading-reloading. 


Reloading [solid] after unloading [dashed] Discontinuous porosity 
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Figure 5.8 Porosity profile under cyclic unloading and reloading at time 
t = 5.8 (or 0.8 after unloading). We have used X = 100, y = 0.2 
and m(t) = 1 ift <5, m(t) = -lif5 <¢t < 5.5, m(t) = 1 if 
5.5 <t <6, m(t) = —lLif6<t<6.5.... Solid part is in compression 
along NCL line, dashed part is reloaded along URL line and dotted 
part corresponds to swelling along URL (unloading). Discontinuous 
porosity occurs at the interface of newly loaded region (solid) near the 
top and reloaded region (dashed) in the middle. 

To understand how the phenomenon of the discontinuous porosity occurs, we refer 
to Fig. 5.9 to aid our discussion. For some short time t after reloading it is possible 
that the effective pressure p < p*, where p* is the maximum value in the time-history. 
The new sediment added to the system will go along the NCL line, while the older 
previously unloading sediment will be reloaded along the URL line. The physical 


condition at the interface between the new and older sediments is the continuity of 
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the effective pressure. From Fig. 5.9, we see that the same value of effective pressure 
corresponds to two different porosity values if p < p*. Thus a discontinuity of porosity 
will appear at the interface. If the loading proceeds to p > p*, then this discontinuity 
will disappear. 

The velocity of the interface of discontinuous porosity can be obtained by using 


the jump condition from the weak formulation of (5.8) 


[l= ¢)u*] 


where u*® = —Vk[ 2 + (1— ¢)J. 
P 4 
oi | 
pb OP inp 


Figure 5.9 Sketch map of effective pressure versus porosity. NCL is 
normal consolidation line and URL is unloading-reloading line. 
From the equation ¢u! + (1 — ¢)u’ = 0, the condition of no fluid stored at the 
interface implies the Darcy flow is continuous, which is equivalent to the statement 


u® should be continuous, i.e., [u*®] = 0. Then the above equation becomes 


[0] [0] 


This means that the interface will ‘fiz’ on the solid matrix. But when reloading 


proceeds to p > p* then [¢] = 0, the interface will move off the solid matrix and start 
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travelling down if effective pressure in the previously unloaded sediment below this 


interface has not reached its previous maximum in time-history. 


5.4 Summary 


The nonlinear compaction behaviour of unloading and loading on a basin scale has 
been modelled as a two-branch nonlinear diffusion equation with a switch condition 
(5.23) and (5.24). The constitutive relations are extended in terms of a dimensionless 
parameter y which is a slope ratio of URL line to NCL line. 

In the case of slow compaction (A << 1), the behaviour is relatively simple, and 
the loading/unloading at the top has very little influence on the porosity profile. The 
behaviour of large \ case (fast compaction) is more complicated and of more interest. 
A downward travelling interface is generated whenever a switch occur between URL 
and NCL branches. The velocity of the travelling interface depends on the slope ratio 
y and decrease with time t. In the lower region, the porosity profile is essentially the 
same as that of constantly increasing loading. 

In the case when newly loaded sediments adds at the top of unloaded sediments, 
a discontinuity of porosity may occur for a very short time. The new sediment 
added to the system will go along the NCL line, while the older previously unloading 
sediment will be reloaded along the URL line. The physical condition at the interface 
between the new and older sediments is the continuity of the effective pressure, which 
corresponds to two different porosity values when p < p*. Thus a discontinuity of 
porosity may occur at the interface. If the loading proceeds to p > p*, then this 


discontinuity will disappear. 


Chapter 6 


Diagenesis: First Order Model 


In the previous chapters, we have mainly investigated the porosity evolution due to 
mechanical compaction. In this chapter, we will analyse the effects of diagenesis on 
the porosity evolution, and show how the model suggests radically different styles of 


behaviour in the distinct limits of slow (A << 1) and fast (A >> 1) compaction. 


6.1 Simplified model equations 


It is clearly seen that Rk, always appears as a combination in the above model 
equations (2.51)-(2.57). It can be easily rewritten as 


ee! 


Rk = exp|G(O > O.)| and 0. = BOR’ 


(6.1) 


where the new parameter 0., which replaces R, is a dimensionless critical temperature 
(with reference to the surface temperature). In the following discussions, we will see 
that the diagenesis reaction virtually takes place in a region called the diagenetic 
window, at a depth of ~ ©,., with its thickness controlled by (. 

From the typical values of model parameters (Smith 1971, Eberl & Hower 1976, 
Lerche 1990, Audet & Fowler 1992) FE, = 60 kJ/mol, T, = 90°C, Ty = 300 K, we 
have G + 2.3,0, % 2, X = 1 and R & 0.01 for d ~ 1km. An initial porosity of 
¢o = 0.5 for pore water at the top of the basin is used by other authors (Smith 1971, 
Sharp 1976, Bethke & Corbet 1988, Audet & Fowler 1992). Initial porosities 0.2 for 


montmorillonite, 0 for illite and 0.3 for quartz are used in our computations. 
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For the convenience of discussing the main effects of diagenesis, we can simply take 


AH = 0, a = 0, and } = 0 in these equations without loss of generality. Based on 


the work of Smith (1971), Sharp (1976) and Audet & Fowler (1992), we adopt the 


following constitutive functions 


P = In(0/@) — (¢0 — 9), 
k = (/0)”, m = 8, 


K =(K,/K,)*”, K,/K, =0.3, 


mh, = 1, 


(6.2) 
(6.3) 
(6.4) 


(6.5) 


We notice that @,. can be determined if we know ¢;, dm, ¢ since ¢.+¢;+¢mt+¢ = 1. 


Thus the quartz equation can be eliminated. Inserting these constitutive relations, 


letting 6 = ¢), and using Darcy’s law (2.56) and the force balance equation (2.58) to 


obtain u’,u!. we finally have the simplified model equations 


Equations for volume fractions 


Ft = e994, — 2 Cibm(l — BRS — Uh 
—(5 — 1a, em), 
Be = (1 = ay)e"- Pn — fell — LESS — 1} 
~(8— 1)a, 02). 
ae = Age (hdl - oS - 1) 
+a; 5e%O-9)g,,, — (6 — 1a APE al 
Temperature equation 
[a(t = 4) FITS = ARE) — a1(6 — a) Gn 
(5 — aww — (a AKC - 8)*15 5S — IS 


(6.6) 


(6.8) 


(6.9) 
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with boundary conditions 


Od = 001 7 
and 
0; na ~j0; Oo a 0, 

: 2 10¢ 
Be TDN a re eee at z= h. (6.11) 

z& 
ie | * ehO-82) 4 de, (6.12) 

0 


where $; = $i, dm; %; and Geo + Gio + Pmo+ho = 1. 
It is based on these equations that the moving boundary problems will be solved 
numerically by using the predictor/corrector implicit finite-difference method pre- 


sented by Meek & Norbury (1982). 


6.2 Diagenesis with slow compaction \ << 1 


From the generalized mathematical model, we notice that the temperature equation 
(6.9) is only weakly coupled with porosity via heat conductivity. If the heat change 
during diagenesis is negligible, then temperature essentially evolves in rather an in- 
dependent way. Thus if A >> 1 the temperature distribution can be treated as a 
prescribed function. In order to investigate the main features of diagenesis, it is con- 
venient to first study the case of \ << 1 (slow compaction) with nearly steady-state 
temperature distribution (A >> 1). When A >> 1 so that conduction is dominate, 


then (6.9) becomes approximately 


ome) 
—S = Ll 
72 0, (6.13) 
so © can be written as 
Ld ae (6.14) 
K 


For the case of A << 1, temperature increases mainly in the boundary layer near the 
basement and is normally not high enough to switch on the diagenetic reaction. The 


temperature distribution is usually very close to steady state in most geological cases 
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of interest. Therefore, we will mainly concentrate on the case of a linear temperature 
distribution. 

The numerical results are given in Fig. 6.1 with a, = 0.15,0, = 2,t = 5. It 
is clearly shown that ¢,, decreases very rapidly in a region of temperature near the 
critical value ©,. The excess pore, overburden and hydrostatic pressures are also 


given in Fig. 6.1. 


1 
illite mont water 
0.8 
N 
£06 be 
o Effective pressure <0 
A 
3 
< 0.4 
oO 
Yn 
0.2 overburden 
hydrostatic 
0 
0 0.2 0.4 0 0.5 1 1.5 2 
porosity scaled pressure 


Figure 6.1 Porosity and pressure profiles with diagenesis (A = 0.01). Z 
is scaled height. The porosities of montmorillonite and illite are marked 
with ‘mont’ and ‘illite’, respectively. The negative effective pressure in 
the diagenetic region is physically unacceptable. Hydraulic failure will 


occur to keep effective pressure non-negative. 


From the numerical results, we understand that @ © ¢0, dz & 0, i.e, u® & 0. This 
means the terms 0(¢,,u*)/Oz in the @» equation, O(¢;u*)/Oz in the ¢; equation are 
negligible. The temperature distribution is approximately (with K= 1) 


© = h(t) — z. (6.15) 


The dimensionless parameter a; represents the effect of the water content released 
during montmorillite diagenesis (typically a, = 0.1). Therefore, it is reasonable to 
assume a, << 1 in the following analysis. The fact that @  ¢p and w << 1 in equa- 
tion (6.11) for the case of small A suggests that h = 1. With these approximations, 


the model equations can be written as 


Obm — —Alh(t)—2-9c] 
BE Ones 6.16 
; é€ ( ) 
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06 _ |, A[h(t)—z-9.] 
- | 2-Odld., 1 
mt aN paz Tue @ (6.18) 


where A’ = A(1 — ¢9)?/¢9. The boundary conditions are 


Q i 0; Om — mo, OF a 0, at z= h, (6.19) 
dz: =, at z=0. (6.20) 
The solutions for equations (6.16), (6.17) and (6.18) can be easily obtained in the 
approximate form 
1 
dm = Pmoexp{— Bern 92} (6.21) 
6: = (1 — a1)(Gmo — bm); (6.22) 


ob = bo — bo V4Nt lerfe 
8a,h oe 


1 p22 2 krz ph kr€ 
EL a i ec — f BIR(t)-E-@clogg @™8 ge 6.23 
“yi D 721 e kcos oh Jo e cos dé (6.23) 


+ 


Solutions (6.21) and (6.22) will be compared with the numerical results later in Fig. 
6.2. It is worth pointing out that solution (6.23) implies that ¢ > ¢o in a narrow 
region near z = h—©,. This is physically unrealsitic which consequently results 
in an interesting phenonemenon known as hydraulic failure, which occurs from the 


diagenetic region up to the basin top (see Fig. 6.2). 


Hydraulic Failure 

From Fig.6.1, we notice that the porosity @ in the diagenetic region can exceed its 
initial value @o; this is physically unacceptable since the effective pressure p. < 0, but 
in reality p. should always be nonnegative. In fact, if p. becomes negative, we expect 
that hydro-fracturing will occur to keep the effective pressure nonnegative. If we 
impose a condition p > 0, the numerical results will ensure that @ < ¢9. But then the 
permeability & will not take the form (¢/¢9)”, and should be determined in another 
way. Hydraulic failure will behave in such a way that an increased permeability 
Kae Will make the fluid drainage balance the water generation to satisfy the physical 


condition p > 0. 
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Since @ © ¢o, o: & O in the fractured region, we have 
w= —(1 _ 0) AK frac: 


Mass conservation implies that 


Okfrac 


= a,behlr)-2-Selg 
De aoe Q 


(1 — bo)? 
Integrating this equation from 0 to z, we have 


a0 i Bln(t)-2-6c] 
rac — 1 : Hee : 
ke = (1 = $)2A 5 e€ o) z 


By using the solution for ¢,,, we have 
Kfrac © 1 for z< h(t) — ©, 


and 


a10¢mo 
(1 — $o)?A 


This means that hydraulic failure can only occur when t > t, & O,¢. 


it 
Rtrac & 1 + exp{— Relies al} for z>h(t)—©,. 


lambda=0.01, Thetac=3, n=2, t=5 Hydaulic failure 
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Figure 6.2 Comparison of the analytical solutions (6.21) and (6.22) with 
numerical results. The permeability kgac (6.28) (dotted) resulting from 
the hydraulic failure increases rapidly at the diagenetic region where 
porosity changes dramatically. The hydraulic failure can occur from 


the diagenetic region up to the basin top. 
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(6.24) 


(6.25) 


(6.26) 


(6.27) 


(6.28) 


The comparison of the analytical solutions (6.21) and (6.22) with numerical results 


is shown in Fig. 6.2 (left figure). The permeability kg,a- (6.28) is also shown in Fig. 
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6.2 (right figure). The permeability ka. (6.28) (dotted) resulting from the hydraulic 
failure increases rapidly at the diagenetic region where porosity changes dramatically. 
It is clearly seen that the hydro-fracturing develops mainly from the diagenesis region 


to the top surface. 


Changing O, and 3 


From the numerical results and the above solutions, we see that the diagenetic 
window is essentially controlled by the depth parameter 0, and the shape parameter 
G3. Fig. 6.3 shows the effect of changing these two parameters. The numerical and the 
analytical solutions are virtually the same, thus we only show the analytical solutions 


in Fig. 6.3. 


Theta_{c}=1.5, 2.5, 3.5 [beta=2.3] beta=2, 5, 8 [Theta_{c}=3] 


illite 


™... Theta_{c}=1.5 


0 005 O1 015 0.2 0 005 O1 015 0.2 
porosity porosity 


Figure 6.3 Changing 0, and @ in the case of t = 5 for \ = 0.01. The left 
figure shows the effect of changing 0, = 1.5, 2.5, 3.5 which only shifts 
the position of the diagenetic region and does not change the shape 
of the porosity profile, while the right figure shows that the change 
of 6 = 2, 5, 8 will dramatically change the thickness of the diagenetic 
region with its central position (@,) fixed. 
It is clearly seen that the change of ©, does not change the shape of the diagenetic 


region but does change its position, while the change of (3 only affects its shape. 
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6.3. Diagenesis with fast compaction \ >> 1 


In this case, the numerical results are shown in Fig. 6.4 with the values of R = 0.01, 
a, = 0.1, t= 5 for a linear temperature distribution O = h — z. It is worth pointing 
out that the curves of smectite (or montmorillonite) and illite in Fig. 6.4 correspond 
to their normalized volume fractions in the solid, namely, solid volume fractions, 
which remove the effect of the change of porosity @ due to compaction. The solid 


volume fraction is related to the real volume fraction by 


1— ¢o T= 


Pm = mT g or bm = PmT— a (6.29) 
Lig DG = L=2@ 
Rad 57 8 b= Pn (6.30) 


illite , 
Dotted: No diagen. 
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Figure 6.4 Porosity profile with diagenesis (A = 100). Z is scaled 
height. The solid fractions (6.29) (smectite or montmorillonite) and 
(6.30) (illite) are used in the left figure. The dotted curves correspond to 
the case of no diagenesis or a, = 0. We see that mechanical compaction 
is the most important factor controlling the porosity evolution and 
diagenesis is only of secondary importance. 
We see that the top region is always nearly in equilibrium, only the middle and 
lower regions are dynamic. It is clearly seen that the diagenesis reaction is essentially 


taking place in a small region, the diagenetic window, which is located in the region 
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where the temperature is nearly at a critical temperature. This figure presents a 
more clear and full view of the compaction evolution during diagenesis. Pore water 
pressure is enhanced by the water released during diagenesis. From this figure, we 
understand that the mechanical compaction is the most important factor controlling 
the porosity evolution, the diagenesis process is also a very important factor, but it 


is in the secondary position. 


Analysis 


From the governing equations, we see that the ¢,,, and ¢ equations are closely coupled. 
Once the solutions for these two equations are obtained, then the solutions for ¢; and 
¢ can be easily determined. For the convenience of analysis, we can take IX (o) = 1. If 
A >> 1, the temperature distribution can approximately be treated as the prescribed 
function 6 = h(t)— z. From the numerical results, we understand that the diagenesis 
reaction is taking place in a narrow region below which the reaction is fully completed, 
and above which the reaction has not switched on. Rewriting the ¢,, and ¢ equations 


(6.2) and (6.3), we have 


Obm _ _ psiht)-2-Edg _ 5 9 7 Pymy (y — gy LOG _ 
z Blh(t)—2z—O¢] 
~(6 — Nae Jo on uae (6.31) 
ODD epOicehs - elon 
= z A[h(t)—2z—-Oc] 
+a, deFlh)-2-8dlg = (6 = 1)ay O{(1 ¢) Jo Smt wm (6.32) 
with boundary conditions 
O 
a -$=0, ab z=, (6.33) 
Om = Pmo, Q —= Qo; 
= Pym _ gy Le _ _ * 4 eblhit)-2-04] 7 
h=1+ AC (1 Ne a 1] + (6 ar [ Ome dz at 2 =f. (6:34) 


From the numerical results and the above governing equation, we understand that 
there exists a transition defined by ¢* = ¢9 exp|—In A/m]. In general, we have ¢* ~ 


O(1), which implies that m >> 1. As the volume fractions ¢,, @; change dramatically 
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in a relatively narrow region, we will assume that G >> 1, and in fact 6 ~ m (or 
B/m = A with A= O(1)). Noticing that the typical value of a; is about 0.1, we can 
also assume that a; ~ O(+) << 1, which will make ma; = O(1). 

From the asymptotic analysis in the case without diagenesis (in Chapter 4) and 
the numerical results, we can easily find that there exist two critical times t* (defined 
as before) and t, corresponding to two typical basin thicknesses h* = II and h, = ©,. 
Normally, t* < t.. For a short time t < t*, the porosity decreases nearly exponentially 
with depth, which means that the compaction is essentially at equilibrium. As time 
increases to a critical time t*, compaction becomes non-equilibrium although the 
diagenetic reaction has not been switched on. As the process proceeds to another 
critical time t., diagenesis comes into play, then we will naturally expect that the 


behaviours may be different in these different cases. 


Short time behaviour (t < t*) 


For a thin layer or short time, we have h — z < II < ©, and exp|G(h — z — 9,)| << 1 


when @ >> 1. Equations (6.31) and (6.32) become approximately 


00m Ce ON 10¢ 
SAR ay oe dm(l — PNB: 1s (6.35) 
and 
ey rat) (1 — 9) Kon (6.36) 


This second equation is exactly the same equation as equation (4.1) we solved in 
Chapter 4. As A(@/¢9)'™ >> 1 still holds, we still can get the Athy-type solution (for 
the leading order) by following the same perturbation procedure as discussed earlier 
in Chapter 4. For simplicity and clarity, we will only repeat some parts of the analysis 


to refresh our solution procedure. The solution for equation (6.36) is 
& = do exp|—(h — z)], (6.37) 


As time increases, porosity @ decreases, but the dramatic decrease of (¢/¢0)™ 
if m >> 1 will cause the perturbation expansions only to be valid if Ak >> 1 or 


b < o* = doexp[—+Ind]. In addition, exp[G(h — z — 9,)] << 1 if @ is relatively large. 
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Now if ¢ > ¢*, (¢/¢*)™ is exponentially large, and exp[G(h — z — ©,)] is exponentially 


small, therefore 
(ey Ga Ooo =m) (6.38) 
using the boundary condition at z = h. We still have 

@ © doexp[—(h — 2)], (6.39) 


whence ¢, © —hd,, and an improved approximation (4.65) to the @ equation (4.1) 


therefore becomes 


(Ge) (1-4) hoe — 1) © h(¢o — 6) — (1 — $0)(1 — A). (6.40) 
To obtain the solution for ¢,,, we change variable by defining 
a dm(1 >. do) 
® = ei a (6.41) 


which is the solid fraction of montmorillonite (or smectite). Combining this with 


equation (6.36), (6.35) becomes 


O® Oo Og @ mr l d® 
(= #5 - OF = OF - (SIs — (6.42) 
By using equation (6.40), we have approximately 
h(1 — @) — (1— 0) 
&, (=) ia (6.43) 
with a boundary condition 
D = dbo. (6.44) 


The characteristics of equation (6.43) imply that 
&’=0, or =x. (6.45) 


In order to obtain the solution for ¢;, we add the ¢; and ¢,, equations to eliminate 
the source terms, so that we have 


05, (du) 


a a, =O With b= b+ (1—a)dm. (6.46) 


The same procedure applies to @ by changing the variable 6 = (1 — ¢o)/(1 — 4) 
leading to ® = (1 — a1) bmo- 
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Transition and solution below transition layer (t* < t < t.) 


The above approximation, however, becomes invalid in the transition region when 


h — z = II or below the transition region and specifically we define 


2 hall a + go 
m m 
6 = dexp[—(-Inm + W)], (6.47) 


whence it follows by a matching principle that V ~ € as € — oo. W satisfies the 


equation 
(—hUe + —W1)b.0exp[—(W — Voo)] 
_— Op ayy pg 2 ee 
= pele" (1 — dooexpl =U — Woo) (Le — DI; (6.48) 
or 
“hbo = wera — b:0)"(Ye = 1]; (6.49) 
where we define 
boo = é*exp[_(—Inm + Uqo)], (6.50) 


and W,, will be defined below. As discussed before in Chapter 4, we still have 


1-0 (1- $00)” vn 


— 51 
h as = (6.51) 
In terms of z, the equation for V, (6.48) is then 
Wedacexpl(W— Wae)] = A fe"{1 ~ dacexpl—(W— Vac). — DV] (6.52) 
(PooexP| ol] = aele exp a ate ; 


Using m >> 1 and exp|(W — U.,)/m] = O(1), the above equation becomes at leading 
order, 


oor + (1 — $o0)e*U, = 0, (6.53) 


The initial data for (6.53) is 
W=W,(7T) when z=0,t=7, (6.54) 
where, if h = 11+ 4Inm at t = to (& #*), then 


Vy(to) = 0, (6.55) 
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and we choose W;,(7) in order that UV, = m at z = 0. The solution is easily found to 
be 
Meine 


vV=lIn 
Tp mltaE ty 


). (6.56) 


This satisfies the boundary condition on z = 0, moreover, we see that (07/0z7)e” = 0, 
so that the diffusion term in (6.52) is identically zero. Therefore (6.56) should give a 
uniform solution to O(1/m) for V in z <h—TII. By matching, we still have 
oo(h — I) 
(1 — Go0)?(t — to) 


All these solutions are essentially the same as those we obtained in Chapter 4. It is 


Uy =In| ] + 0(1). (6.57) 


worth pointing out that the solutions obtained so far are only valid for t < ¢. and 
there is no reaction involved. If the reaction comes into action, then the boundary 
condition (6.54) is no longer valid because the base z = 0 is not reachable. The 


boundary condition will be modified accordingly in the following subsection. 


Intermediate region (t > t., h-O.< z<h-—TID) 


If diagenesis is taken into account, the improved approximation to the solution (6.40) 


of the ¢ equation should be modified to include the diagenesis term. It becomes 


a 
o 
Now the term exp|[3(h — z — 9,)] is still small in the region z > h(t) — ©, we still can 


()"(1 = 6)'(a— = 1) & h(a -— 4) — (1 = Go)[L =k = a1(5 — 1) moh] (6.58) 


expect there exists a similar transition region (z ~ h(t) — II), and we will have the 
same equation as (6.53) below the transition region, but now the boundary condition 


is different because the base is not reachable. The boundary conditions are 
Vowvy, as z—-h-—Il, (6.59) 


UV, a z>-h-OQ, (6.60) 


where ¢. = ¢*exp[4+(—Inm + W,)] and $x = ¢*exp[+(—Inm + V,.)] are to be deter- 


mined later. The characteristics of equation (6.53) are 


b =0, p= otek (6.61) 
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whence 
_ _ (1 a ne Ww 
=v... 2= ae a ({—7)+A(r) — ©,, (6.62) 
which becomes, by using the boundary condition (6.59), 
ie 
h-W= cae (t—7)+h(r) — ©. (6.63) 
Eliminating 7, we have 
= Tk=s\p., 
U =U, lt - (h 210 ). (6.64) 


(I= Oye )*(ere ere) 
This solution will determine ¢, as z — h — II when ¢, (or V,) is known. 
In order to determine h, we rewrite equation (6.58) in terms of € and t, and match 


it to equation (6.48); we thus find that 


ies 1— dp Meg (1 = $0)" ow. nag (6.65) 


Clearly, if there is no diagenesis R = 0 (O, — oo) or diagenesis without water release 
(a, = 0), the above expression will degenerate into (6.51). 
Reaction region (t > t.) 


In the region z ~ h(t) — ©., the term exp[G(h — z — 9,)] will not be small, we can 


expect that there will exist another transition in this reaction region. We define 


fi =o BG 
z=h-O, B + 8 
ie *exp[—(—Inm +). (6.66) 


By changing variables in this way, we have exp|[G(h — z — 9.)| ~ @ and we thus 
balance the terms in the governing equations (6.31) and (6.32). By using the chain 
rules 0, = GOc, O; = Oj — Bhd, WV and @¢,, satisfy the equations 


1 , (1-45)? 8, 6/8 may © \o~6 
gu =hVe= BoC [e"(— We =1)|+ me” —(6-1)(1—¢h)leSmn;, (6.67) 
10m ;8bm — (1— 4%) O(bme®) 2 a Sues 
pee 5 Ome — (6-1) zelOm fone dc], (6.68) 
where 
ie é*exp[—(—Inm +0), (6.69) 
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and WS, will be determined later. 
If we use the conditions m ~ 3 >> 1 and ma, ~ 1, then the above equations can 


be written approximately as 


B 


(age) Oye pO ms 
—hWe — — e€ We —1)] = —[6 — (6-1) — 6% )le Co abm, 6.70 
— SSE Sle Suc — = SE (6-0 SNe Som (6:70) 
- Obm -¢ 0 S = 
Th = Om = =] B-lPm m ¢ ’ Gk 
ape = ome — (5— Van zelbm J dme Sd) (6.71) 
where we see that the approximation of G ~ m and 2 ~ O(=) are appropriate. 
The far field matching conditions are 
Loe 
dm > bmol ), Ur, as (> ~, (6.72) 
1— % 
VWs as € > —ov, (6.73) 


where the factor (1 — $S,)/(1 — @o) is due to the effect of porosity change from ¢ to 


bf. Since ay << 1, [S., dme~Sdz = O(1), then the ¢m equation becomes 


:Obm 2 ¢ 
hae = bme$. (6.74) 


Integrating this equation and using the matching condition, we have 


dm = Bboy —op2 )expl— 5 €~4. (6.75) 
Substituting this solution into the V equation and integrating from —oo to ¢, we have 
iw ol oF y.1))-B 2 me [5—(5-1)(1-6%)} no — Jexp[— Fer 
(6.76) 

where 
B= [hve - (= $50)" (vs) (6.77) 


The solution for V can be written in a quadrature. Clearly, when ¢ — -+too, we 


obtain the jump condition in the diagenetic region 


[—hw fs (1 aoe) om, 
= hmay 


MAS — (6-1) - 65) bol 


), (6.78) 
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which determines WV, in term of US. This gives a shift in the porosity @ outside the 
reaction region. If we define the thickness of the reaction region as the distance that 
dm Change from 90% to 10% of the initial value, then the dimensionless thickness of 
the reaction is about 


6.3, we have In[10G]/G ~ 1.5, 0.8, 0.5 for G = 2, 5, 8. 


me which is clearly seen in the numerical results. In figure 


Solution below the reaction region 


To obtain the solution for VS, as ¢ — —oo, we write the equation for V in terms 
of z. By using the solution (6.75), we see that the source term due to diagenesis is 
virtually negligible in the region z = 0 to h — O,. Then the equation for V becomes 
(at leading order) 

OU + (1 — 6 )?e" Ws = 0. (6.79) 


Co 


Following the same procedure as before with the boundary condition V, = m at the 
base z = 0, then we still obtain (6.56), and finally we have 

b5(h = Q,) 
(Page) eat) 


which completes the solution procedure. 


we, =In ]+0(1), (6.80) 


Summary and Comparison 


The solution of equation (6.32) with boundary conditions (6.33) and (6.34) consists 
of a near equilibrium solution (6.39) in the upper region, a transition given by (6.49), 
an intermediate region (6.64), a reaction region (6.78) and the solution below the 
reaction region (6.80). 

Solution (6.80) gives US (and @& through (6.66)), (6.78) determines WV, (and ¢,) 
in terms of VS, (6.64) gives VU. as z — h—II, and (6.65) provides an equation which 
determines the evolution of h(t). 

The comparison of the solutions (dashed) with the numerical results (solid) is 
shown in Fig. 6.5, Fig. 6.6 and Fig. 6.7. It is clearly seen in figure 6.6 that the 


agreement gets better as @ becomes larger. 
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Comparison with solutions 


Dashed: Sol 
Solid: Numeric 


0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 
Porosity 


Figure 6.5 Comparison of solutions with numerical results. The values 
of A = 100, a4, = 0.1,¢ = 5,8 = 2.3,0, = 3 are used. The dashed 
curves are calculated from solutions (6.40) (top), (6.76)(middle) and 


(6.80) (lower). 


Changing beta [=5, 8] 
0.5 1 


Dashed: Sol 


0.45+ Solid: Numeric 


N 0.357 


0.3/7 


15 ‘ 02 0.25 

orosity 
Figure 6.6 Comparison of solutions with numerical results. Parameters 
as for figure 6.5, but for different values of @ = 5, 8. The dashed curves 


are calculated from solutions (6.76)(middle) and (6.80) (lower). 
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mont 


upper solution 


reaction region 


Dashed: Solution 
Solid: Numeric 


illite 


0 0.1 02 San 03 0.4 0.5 
Figure 6.7 Comparison of solutions with numerical results (A = 100, t = 
5). The solid fractions (6.29) (smectite or montmorillonite) and 
(6.30) (illite) are used in this figure. The dashed curves are calculated 


from solutions (6.45) (top) and (6.75) (lower). 


6.4 Application 


The data analysis given by Abercrombie, Hutcheon, Bloch & Caritat (1994) from 
oceanic and sedimentary basins shows that burial history has significant influence on 
the Smectite-illite (S-I) diagenetic reaction. In a slow burial environment, the S-I 
reaction may begin at temperatures as low as ~ 50°C, and reaches completion by 
~ 90°C; while in a rapid burial environment, the S-I reaction may not begin until 
temperatures as high as ~ 120°C, and may not reach completion until ~ 150°C. From 


these results, we understand that 


e Diagenesis takes place at lower temperatures or shallower regions in the fast 


compaction process (A >> 1) than in a slow compaction process (A << 1). 


e The diagenetic process is essentially constrained to a narrow region (a diagenetic 
window ) with a temperature range ~ 30°C or equivalently over a depth range 


of ~ lkm. 
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By using the present model and the solutions obtained so far, we can explain these 
phenomena. From the definitions of the parameters, we find that the depth to the 
centre of the reaction window; d,, is 
7 KoRTy Ms 


UO? ares. Id 


(6.81) 


This clearly means that the higher the sedimentation rate, the higher the critical 
temperature of diagenesis, the deeper the diagenetic region, and vice versa. A change 
of 2 orders in sedimentation rate will cause a shift of d. by 2 (equivalently ~ 60° 
C) (with other parameters unchanged). In addition, the thickness of the diagenesis 
region dg; is the order of m0) d, A typical value of @ % 2.3 gives dgy © 1.36km 


(with d = 1 km), or equivalently a temperature range of ~ 40°C. 


Chapter 7 


Diagenesis: Dissolution and 


Precipitation Model 


The smectite-to-illite transformation is the most important process during shale di- 
agenesis. The mathematical model presented in the last chapter is a first-order dia- 
genetic reaction (dehydration) model in which the geochemical compositions of pore 
fluid are not taken into account. The main factor included in the model is tempera- 
ture. In reality, diagenesis is far more complicated and takes place via the dissolution 
of smectite in pore water and the subsequent precipitation of illite involving the inter- 
actions of many mineral species. This chapter’s purpose is thus dedicated to extend 
the first-order dehydration model in Chapter 6 to a more realistic reaction-transport 


mathematical model with a more detailed analysis in some practical cases. 


7.1 Introduction 


Diagenesis is observed world-wide in sedimentary basins. The close spatial and tem- 
poral correlations between smectite disappearance and illite formation imply the ex- 
istance of the smectite-illite reaction. Such a smectite-illite (S-I) reaction is one of 
the fundamental mechanisms in clastic diagenesis. The reaction has received much 
attention but the nature of both the illite/smectite (I/S) mixed-layer and the reaction 


mechanism are still under discussion, and many experiments have been carried out 


100 
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to investigate the kinetic features of the S-I conversion (Eberl & Hower 1976; Bethke 
& Altaner, 1986; Huang, Longo & Pevear, 1993; Abercromie, Hutcheon, Bloch & 
Caritat, 1994). 

Many authors write the overall S-I transformation with K-feldspar as the following 


reaction 
smectite + K-feldspar — illite + quartz (aq) + interlayer water. (7.1) 
One detailed example of this symbolic reaction is 
K AlSiz0g(K-feldspar) + 2h9.3Al1.9(Sis)O19(OH)2 - 4.520 (K-smectite) 


2K 8Ali.9(Alo5Si3.5)O19(OH )s (illite) + 9H2O (interlayer water) + 4.57%02(aq), 
(7.2) 
as given by Abercromie, Hutcheon, Bloch & Caritat (1994). 
Recently, Huang, Longo & Pevear (1993) systematically analysed experimental 
and field data and derived the conversion rate 


2 = KDR S? or T= HTN 2), (7.3) 


where S, J are relative fractions of smectite and illite in the 1/S mixed-layer ([+S = 1), 
[A*] is K* concentration in the fluid, k(T’) is the reaction constant depending on 
temperature 7’. 

We see that the experimentally derived reaction rate (7.3) is second-order with 
respect to smectite and first-order with respect to A* concentration. This empirical 
relation can be easily obtained via the law of mass action for the reaction (7.2). 

Potassium cation concentration has an important effect on the reaction rate. At is 
mainly supplied by the dissolution of K-feldspar. The characterization of K-feldspar 
dissolution rate is essential for the accurate description of the overall S-I process. Ex- 
perimental investigations of feldspar dissolution rates have been performed (Busen- 
berg & Clemency, 1976; Helgeson, Murphy & AaGaard, 1984; Chou & Wollast, 1985; 
Hellmann, 1994; Gautier, Oelkers & Schott, 1994). Although the dissolution process 


of K-feldspar in natural environment is nearly at equilibrium with a temperature 
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range of 0 ~ 150°C, 
K AlSiz0g + 2H2O = Kt 4+ 3S9i02(aq) + AI(OH);, (7.4) 


most of the experiments have been carried out at far from equilibrium conditions with 
a temperature range of 200 ~ 400°C. 

Helgeson (1968) and his coworkers (Helgeson, Garrels & MacKenzie, 1969) devel- 
oped the first model to consider water-rock interaction as a system of coupled dis- 
solution and precipitation. In their model, dissolution reactions of primary minerals 
(smectite, K-feldspar) are treated as irreversible processes, while partial equilibrium 
with respect to the secondary phases (illite, quartz) is assumed. As pointed out by 
Helgeson (1979) and Steefel & Cappellen (1990), the assumption of partial equilib- 
rium is only justified where the rate of precipitation of a secondary phase is faster 
then the rate of dissolution. However, the precipitation of the stable insoluble min- 
erals may be slow even on geological time scales. Therefore, the partial equilibrium 


may be a good approximation to natural water-rock systems. 


7.2 Mechanisms of S-I Reaction 


Extensive studies on clay diagenesis with increasing depth of burial reveals that the 
most systematic evolution consists of the progressive illitization of smectite minerals 


(Chamley, 1989). 


e Firstly, such modifications usually occur at depths exceeding 2 km. In the series 
marked by normal geothermal gradients of about 30°C/km, the process develops 
between 2.5 and 3.5 km, and does not progress beyond a depth of 5 or 6 km. 
This suggests that highly expandable smectite-rich minerals change to slightly 
expandable illite-rich ones over a relatively narrow temperature interval (the 


diagenetic window); 


e Secondly, contrary to what is often believed, diagenetic processes do not notice- 


ably depend on the absolute age of burial series. More important than geological 
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age are certainly the geothermal gradient and the residence time at a diageneti- 


cally active temperature; 


Thirdly, the illitization process does significantly involve A+ which is supplied 
by dissolved K-feldspar. 


In addition, extensive investigations also suggest that some cations like At, Alt, 


needed for the evolution of smectite to illite, are provided through short-distance 
transportation not by long-distance transport processes. The dominance of such 
very short exchanges in shaly sediments is confirmed by high-resolution trans- 
mission electron microscopic observations (Ahn & Peacor, 1986). This means 
that diagenetic sediments are not significantly affected by pore-water migration, 


at least not since their initial compaction. 


e Finally, the argillaceous deposits (in the Gulf Coast) behave essentially as a 
nearly closed system, pore fluid being present in small amounts compared to the 
solid materials, and acting possibly as a catalyst for short-distance ion transport 
and for local clay reconstruction at the reaction interfaces (Ahn & Peacor, 1986; 


Chamley, 1989), which means that the system is nearly at equilibrium. 


Several mechanisms have been put forward to explain the S-I reaction process. 
Two main ones are transformation and dissolution-precipitation. The former mecha- 
nism suggests that the S-I reaction is a transformation process through smectite/illite 
mixed-phase with (a series of) reordering processes of the intermediate mixed-layer 
(Hower et al, 1976). An alternative modification is a solid-state transformation mech- 
anism without mixed-layering. The latter mechanism involves the processes of smec- 
tite dissolution and illite precipitation without mixed-layering. According to high- 
resolution electron microscopic data, the mixed-layering mechanism appears to be 
questionable (Chamley 1989), but Ahn & Peacor (1986) provide a seemingly convinc- 
ing example of a smectite-to-illite transformation rather than a neoformation. The 
first-order dehydration model in the previous chapter is essential a transformation 
model. The fact that diagenesis, which is still imperfectly understood, largely de- 


pends on lithology, fluid pressure, geothermal gradient and pore fluid compositions is 
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the main motivation for us to develop a more realistic reaction-transport dissolution- 


precipitation model in the present work. 


7.3 Model Equations 


The S-I transformation is composed of the following intermediate dissolution-precipitation 


reactions which can be written symbolically as 


Smectite dissolution 


M* (smectite)  [X“] + n[H20], (7.5) 
Illite precipitation 
[At] + [AlOs’—] + f[X*] 3 fr (illite) + [Si0#], (7.6) 
K-feldspar dissolution 
[K-feldspar] 3 [K+”] + [AlOz"—] + s[Si05], (ae 


Quartz dissolution and precipitation 


[SiOZ] 4 [quartz], (7.8) 


ra 

where n, s, f are stoichiometric coefficients and S, L denote solid and liquid phase. [X] 

is an aqueous silica combination in the form such as [—(S%4)O10(OH)2]. [AlOz”—] is 
only a general notation of the combination such as [Al(OH); |. 

Let the molar rates of the above reactions be r1, 2, 73,174 (for forward reaction) and 

r_4 (for backward reaction), respectively. Then, the reaction-transport (by diffusion 


and advection) model can be written as 


o[M] ea 
a + Vw Ml] = —r, (7.9) 
Ae) +V -[u'g[X]] — V - [2DV([X])] = 1 — fra, (7.10) 
all 


OE. +V- [u®[Z]] = fro, (7.11) 


CHAPTER 7. DIAGENESIS: DISSOLUTION AND PRECIPITATION MODEL 105 


AED £7 (wlo[K J -V-[SDVK J =-re brs, (7.12) 

AED 7 fu'o[Ar]] —V-[DV(AM)] =—re tr, (713) 

eS tOe) + V-[u'd[Si02]] — V - [6DV([Si02])] = srs —ratreatre, (7-14) 
a +V-[u'[H,0]] = nny, (7.15) 

Paquarte) yaaa Seka (7.16) 

olf ae ares [us| feldspar]] = —rs, (7.17) 


where [MM], [J], [quartz] are molar concentrations, measured in units of mol m~? of 
rock. [SiOz], [A]... are molar concentrations in units of mol m~® of pore water. ¢ is 
porosity. 

The reaction rates r;,2 = 1, 2,3, 4, —4 are generally complicated nonlinear functions 
of concentrations, satisfying r; = 0 at equilibrium (Dewynne, Fowler & Hagan 1993). 


Their precise form should be determined by experiments. 


7.3.1 Surface controlled or transport controlled 


The kinetics of mineral dissolution and precipitation are strongly controlled by reac- 
tion rates which depend in a complicated way on the solution compositions and surface 
chemistry. A complete formulation of a quantitative dissolution /precipitation rate law 
is more complicated for multicomponent systems. It is usually helpful to identify one 
of the processes, transport or surface attachment, as the rate-limiting step (Berner 
1978, Lasaga 1981, 1984). If the transport process is much slower than the reaction 
rate at the surface of the mineral, then the dissolution and precipitation are referred 
to as transport-controlled, while the opposite case is termed surface-controlled. 

In the case of a transport-controlled process, the surface detachment and attach- 
ment are so rapid that a saturated solution adjacent to the surface is maintained. 
Dissolution and precipitation are then regulated by transport via diffusion and ad- 


vection into the surrounding medium. The reaction rate thus depends on the flow 
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velocity and the degree of stirring (Nielsen, 1964). If transport is by pure diffu- 
sion without advection, then the case is termed diffusion-controlled (Berner 1978). 
In surface-controlled dissolution/precipitation, the surface process is sufficiently slow 
that it can not keep pace with diffusion and advection. The concentration level adja- 
cent to the surface is essentially the same as that in the bulk solution. Flow velocity 
and stirring have a negligible effect on dissolution/precipitation rates. Generally 
speaking, surface-controlled dissolution/precipitation reaction is slower than that by 
transport-controlled process. The two limiting cases are determined by the final value 
of the surface solution compositions. 

Berner (1978) found that most mineral dissolution/precipitation reactions are very 
close to the case of a surface-controlled process. We will see below that the rate laws 
derived from most experimental data are essentially surface-controlled. Therefore, 


the transport effect is included in the model equations but not in the reaction rates. 


7.3.2 Nucleation and crystal growth 


The precipitation process can be described in more detail as nucleation and crystal 
growth. If the concentration is gradually increased, exceeding the solubility with 
respect to a secondary solid phase, the new phase will not form until a certain degree 
of supersaturation has been achieved. Stable nuclei can only be formed after an 
activation energy barrier has been surmounted. Nucleation normally proceeds via 
homogeneous or heterogeneous nucleation. In most cases, however, heterogeneous 
nucleation is the predominant formation process in natural waters since it has a lower 
activation energy barrier than that in the case of homogeneous nucleation. Just as a 
catalyst reduces the activation energy of chemical reaction, foreign solids may catalyze 
the nucleation process by reducing the energy barrier. Phase changes in natural 
aqueous systems are almost always initiated by heterogeneous solid substrates, as 
pointed out by Stumm (1992). 
The free energy of heterogeneous nucleation AG; can be generally written as 


Q 


eq 


AG; = RTIn + AGsurt(7), (7.18) 
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where K,, is the equilibrium constant for the reaction, and @ is the reaction quotient 
for the same reaction. Q/K., represents the solution saturation state (Steefel & 
Cappellen, 1990; Stumm 1992). ¥ is the interfacial free energy. For homogeneous 
nucleation, AG. = YowA. In the case of heterogeneous nucleation, this term should 
be modified because the nucleus is now formed in part contact with the solution and 
in part with the surface of the solid substrate (Cappellen, 1991; Stumm, 1992). More 
generally, we have 


AGsurt = YowAcw + (Fcs — sw) Acs; (7.19) 


where the suffixes CW, CS, SW refer to cluster-water, cluster-substrate and substrate- 
water, respectively. If ¥sw >> Yow, the precipitate tends to form a structurally 
continuous coating on the substrate grain. In this case, the interfacial energy may even 
possibly become negative and the activation barrier vanishes. These considerations 
show that the interfacial energy is of importance in determining the thermodynamics 


and kinetics of the nucleation process. 


7.3.3 Rate laws for dissolution and precipitation 


Most dissolution/precipitation experiments are carried out under far from equilib- 
rium conditions. However, such laboratory data are not directly applicable to field 
observations. Unfortunately, the discrepancies between field estimates and laboratory 
measurements of reaction rates are as large as up to four orders of magnitude. One 
possibility of explaining this difference lies in the fact that not all of the potentially 
available surface in natural systems actually participates in reactions with pore fluids. 
A common implicit assumption in modelling interface-controlled kinetics is that the 
rate is linearly dependent on surface area of which is poorly estimated in spite of its 
vital importance for a better understanding of the reaction mechanism. Coating of 
mineral surfaces by secondary mineral precipitation and associated occlusion of natu- 
ral surfaces may account for the apparent lesser reactivity of natural mineral surfaces 
relative to their laboratory conterparts. However, the extensive etching widely ob- 
served on some silicate-mineral surfaces indicates all portions of the primary mineral 


surface are accessible to pore fluids in spite of secondary precipitation, militating 
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against an extensive surface-covering role for coating (Velbel, 1993). Lasaga (1981, 
1984) and Aagaard & Helgeson (1983) try to bridge the gap and extend the laboratory 
kinetic data into a general rate law that is applicable to natural situations. 

In the case of interest to us, smectite dissolution will normally proceed with respect 
to the nearly amorphous silica solubility 2 x 107? (M) or 120 ppm at 25°C. Quartz 
precipitates with respect to quartz solubility 1x 10~* (M) or 6 ppm. Illite precipitation 


illite 


eq between the upper limit of amorphous silica 


goes with respect to a solubility c 


amorph 
eq 


solubility ¢ and the lower limit of quartz c%"" according to the thermodynamic 
and kinetic constraints and the activity calculations by Aagaard & Helgeson (1982, 
1983). 

According to the earlier works by Rimstidt & Barnes 1980, Lasaga 1981, 1984, 
Ortoleva, Merino & Sen, 1987 and Huang, Longo & Pevear, 1993, we can generally 


write the reaction rates as 


where k;,7 = 1,2,3,4,—4 are rate constants which are functions of temperature 7’. 
A;,... is the specific reactive surface area (m?/m?°) of the mineral (smectite, illite, 
K-feldspar and quartz), f;(a;) is a function of the activities a; of the jth primary 


species in solution, which is usually assumed to be of the form 
fila;) = si Get (7.21) 
j 


where n, is the stoichiometric coefficient. g(AG;) accounts for the important variation 
of the rate with the deviation from equilibrium (AG; = 0). Lasaga et al (1994) write 


this function as the following form 
g(AG;) = (1—exp(AG;/RT)), (dissolution), 


g(AG;) = (exp(AG;/RT)— 1), (precipitation), (7.22) 


where (¢)+ = maz{0,¢}, AG; is the Gibbs Free Energy of the reaction. AG; < 0 
is for undersaturation, while AG; > 0 for supersaturation. Note that this equation 


satisfies g(0) = 0. At constant temperature and surface area, it follows from this 
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equation that the dissolution rate will be essentially constant at far from equilibrium 
condition (AG; << 0). Normally, for a single species dissolution and precipitation, 


AG; = RT In(c/c-,), then we have 
g(AG;) = (1—$), (dissolution), 


g(AG;) =(S—1)4 (precipitation), (7.23) 


where S = c/Ceq is the saturation ratio, and c,, is the concentration at solubility 
equilibrium. The present ability to predict reaction rates as a function of saturation 
state is still limited. A further consideration is the possibility of fully nonlinear rate 


laws. Rate laws with a functional dependence on AG; of the form 
g(AG;) = (1 — exp(AG;/RT))™, (7.24) 


have been applied most commonly to precipitation kinetics (n; = 2). 


The temperature dependence of k; follows the Arrhenius law 
ky = vjeP/FP (i = 1,2,3,4, —4), (7.25) 


where E; is the activation energy, v; is the frequency factor and R is the gas constant. 
A; is a function of the volume fraction of the mineral. For uniformly packed spherical 
particles with an averaged radius 7;, we have A; = 3¢;/7; or A; x @j. 


For convenience in the following discussion, we can rewrite the rate laws as 
r, = ke exp(—E;/RT) fj, (7.26) 


where ker ! has the unit of s~!. r; absorbs all the other terms (noting that A; « ¢;) 
and has the same units as molar concentration. If r; is written in terms of volume 


fractions, then we have 
density 


A 


See ee 2 
molar weight ee) 


where 7; is dimensionless. The term (density)/(molar weight) can be written as 


Pm/Mm,t = 1,2 and psi/Msi,i = 3,4, —4. 


CHAPTER 7. DIAGENESIS: DISSOLUTION AND PRECIPITATION MODEL 110 
7.4 Non-dimensionalization 


Let dm, i, @f,q be volume fractions of smectite, illite, feldspar and quartz, respec- 
tively. cx, cs; and cx are the solubility limits of [X], [SiO2] and [A *],respectively 


By using the relations between molar concentrations and volume fractions 


[v= PrP, = EE lauarts] = SP, [a0] = Fe, [feldspar] = PPE, 
_ PXOx 1 PsiPsi _ PKOK PAP al 
[xX] = Ma” [SiO2] = “Mos” [K] = a [Al] = ea (7.28) 
we can write the governing equations in terms of volume fractions 
a + V+ (u'¢m) = n=, (7.29) 
+V- (alox) — V-[D(6) Vox] = (ri — ra), (7.30) 
PX 
8 V-(u'ds) = fra (731) 
M. 
ven +V-(u'dx)-V- [D(¢) Vox] = (-r2 + r3)—, (7.32) 
PK 
ys (u'gar) — V - [D(¢)Veai] = (—1r2 + Bees (7.33) 
ot PAL 
78s +V-(u'dsi) — V - [D(¢)Vesi] = (sr3 — ra + 7-44 pees: (7.34) 
Psi 
So V (u's) = (ns nar), (7.35) 
Et V: (wy) = 4rd (7.36) 
768 +V-(u'ds) = 1? (7.37) 


where D(@) is a known function of ¢. 
If we scale u’,u! with m,, z with d, t with d/m,, dx with éx, dg; with d5;, dK 
with dn, kf with ROM 


eet — pep § = 1,9,3,4,—-4, (7.38) 


where $x, @s:,¢K,a are the volume fractions corresponding to the solubility limits 


of [X], [SiOg], [K+], [Al*]. As before, T is rescaled as O = (T — To) Ko/qod which is 
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the dimensionless temperature with reference to the surface temperature J. Then 


the governing equations become (without diffusion) 


Obm naan 
“om + V : (u°dm) = —Ryikiry, (7.39) 
O a a a 
nen + V (u'dx) = — [Riki = fRoakoro), (7.40) 
t ox 
Od; : - 
at + V -(u°d;) = fRek2aer2, (7.41) 
0 1 = = 
nn +V-(udx) = —[-a3Rokofo + agR3k3F3], (7.42) 
OK 
Ooai I as cs — 
+ V+ (wear) = =—|[-a3Rekoho + asaRaksrsl, (7.43) 
Ot PAL 
Odsi I te 1 > - > a3 = sy 
+ V : (u ogi) == [sR3k3r3 = Rakara + R_ak_aT_4 + Roker, (7.44) 
Ot si a4 
Ze V - (u'd) = agRi kif 7.45 
at (wd) = agRikiri, (7.45) 
Obq s = faves Le 
a + V-: (u bq) = a7(Rakara = Fak aren, (7.46) 
O a 
er + V i (u’ds) = —agR3k3r3, (7.47) 
where 
eff g 
Ret Se SA (7.48) 
Ms 
7 ; Exqod 
— ei fos 
k; =e with 6; KoRT2 (7.49) 
ee PmMx PmM; PmMx psiMr 
ac = 
pxMn PiMm PKMy, pKMs; 
pKMa NPmM wy psiM, psiMs 
paMx PwM, pyMs, prMsi en) 


We see that this present model is essentially controlled by ten dimensionless 
parameters R; and 3;,4 = 1,2,3,4,-4. a; = O(1). According to Rimstidt & 
Barnes (1980), Aagaard & Helgeson (1982,1983) and Stumm (1992), dx ~ 12 x 
10->(amorphous), dg; ~ 6 x 10~®(quartz). ie. éx, OK, 051,041 << 1. These condi- 
tions suggest that the reactions for the aqueous species are very fast, and thus the 


pseudo-steady state approximation is valid. 
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7.5 'Two-step Case and Dehydration Model 


In order to investigate the main features of the dissolution-precipitation process and 
verify the validity of the first-order dehydration model discussed in chapter 2 and 
chapter 6, we will now simplify the whole model into a two-step case of one-step 
dissolution and one-step precipitation. With these simplifications, we can analyse 
in more detail the effect of transport and reaction rate and find out the controlling 
mechanism of the processes. We will also see that how the dehydration model can be 
derived from this two-step model. 

In the two-step case without potassium/aluminium/silica activities, ky = k_4 = 
ks = 0, 6x = 1 (at equilibrium), f = 1. This corresponds to the following reaction 
mechanism 


Smectite dissolution 


M* (smectite) 4 [X"] + n[H20], (7.51) 
Illite precipitation 
[X4] 3 I (illite). (72) 


It is clearly seen that R;k; always appear as combinations in the model equations. 


This combination can be easily rewritten as 


_ 1 1 
Rik; = exp[G;(0 — O,;)] and O,; = —n—, (V59) 
Oe I 
ile (O)eff 
Ri = fs : a fee (7.54) 
it 


where the new parameter ©,;, which is equivalent to R;, is a dimensionless critical 
temperature (with reference to the surface temperature Oo). 
The above model equations can be rewritten as the following one-dimensional 


model (without transport) 


Ibm 


ar = TW Rabat, (7.55) 

09x = ( = do) [Rikir = Rokof ol, (7.56) 
ot € 

09% = Ro2(1 —= dg) kara, (7.57) 


ot 
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Og = agdRi ki", (7.58) 


with 


Ds 
and 6 = —. 7.59 
M; Pl ( ) 


where we have used the approximation My = M;, Mm = M;+nMy, Pm = pi = Px = 


Ps: Pw = pl, and e = dx << 1. We notice that the first two equations (7.55) and 
(7.56) are decoupled from the other two equations (7.57) and (7.58). The total mass 
conservation implies that ¢@, + ¢; + 6/6 + €dx =1. 

For typical parameters kl! ~ 1 x 10716 s-1 (Swoboda-Colberg & Brever, 1993), 
KOT ~ 0.2 x 10-165"! (Small, 1993), d ~ 1000m, mh, ~ 0.5 x 10-4 ms7!, Ry & 
0.02, R. =~ 0.004. Here we have used Swoboda-Colberg & Brever’s results (1993) 
that dissolution/precipitation rates measured in the field appear to be a factor of ~ 
200—400 slower than that of the same minerals measured in the laboratory. E, ~ 60— 
80k J mol! (dissolution) (Eberl & Hower, 1976; Lasaga 1984), Ez ~ 90—110kJ mol! 
(precipitation) (Small, 1993) (£2 > F,) correspond to (3; ~ 2.3 — 2.8, G2 & 2.9 — 3.5 
and 0.) % 2 < O.o & 2.15. 


7.5.1 Degeneration to the dehydration model 


Now we will show that how the present two-step model can degenerate into a first- 
order dehydration model and thus verify the validity of first-order model we have 
discussed in chapter 2 and chapter 6. 


Since € << 1, the steady-state approximation for [X] can be used. Thus we have 


Rikiry = Rokois ~ 0. (7.60) 


If we assume ky = ko = ky (R2 = Ri = R), and use the following rate functions 
T9 =T) = Oras (7.61) 


then we obtain (with transport) 
Om, A(u°br) 

Ot Oz 

Oo; , A(u*gi) 


at | Oz 


= —Rk,dbm, (7.62) 


= R(1 — ao) krbm, (7.63) 
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do, Ou'd) 
ot Oz 


These equations are exactly the equations (equations (6.2)-(6.4) with ay = a) we 


= ab RKpbm: (7.64) 


discussed before in the one-step dehydration diagenetic model. Therefore, the same 
solution procedures will give the same results as those in Chapter 6. This first-order 
approach is consistent with the previous attempts (Eberl & Hower, 1976; Bethke & 
Altaner, 1986), and can represent adequately complex reactions with a rate-limiting 


reaction step (Lasaga, 1981; Velde & Vasseur, 1992). 


7.5.2 Effect of transport 


From the above subsection, we see that the assumptions (k, = ko, 7 = 71 = @m) are 


very specific. To be more realistic, we use the following rate functions 


ry = dm(ds > @x)4, 9 == (dx = ds)+, (7.65) 


where ¢, = $;(Q) is the solubility (of [X]) which is a known function of temperature 
or equivalently a function of time t and depth z. ¢,(z = A(t)) = 1. @, usually 
increases as © increases. 


The first two model equations with transport are 


O%m O(u?om) 


a ae —Rikidm(os — dx) +, (7.66) 
O O(u! 1 0,~0 1- = - 
Ox 4 ex) 8 (BOOX) = 7" RF ble — x)4 — Rakaldx ~ 0)4) 
(7.67) 
where 
msdTH 
Pe= Do (7.68) 


is the Peclet number. Dp is the diffusion coefficient, and 77, is the tortuosity. For 
the typical values of m, ~ 0.5-"ms~', d ~ 1000 m, Do ~ 1079 m? 871, 74 ~ 8, 
Pew 15. 

Reaction without transport 


In this case, the equations for dm», dx become 


bm 7s —Rikibm(ds + bx )+, (7.69) 
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peel INC ee Pena eC eee (7.70) 


€ 


Since € << 1, the above equations are similar to the model equations for enzyme 
kinetics. Therefore, Michaelis-Menten’s pseudo-steady state approximation applies. 
In other words, the reaction for ¢x is so fast it is more or less in equilibrium at all 


times. Mathematically, we have 


Rikibm(ds — x)+ © Reoko(ox — os), (7.71) 


which holds exactly only if 6, = ¢,. Thus, we can simply look for a perturbation 
ox = ost GY +.... (7.72) 


Subsistituting into equations (7.69) and (7.70), we have 


bm = —Rikiebm(—oY) 4, (7.73) 
bs = (1 — a0) [Rikidm(—bY)+ — Reka(oY)4. (7.74) 


We can easily see that if oY <0, then 


bp t —— 6, = 0. 7.75 
Pm + 7 ae (7.75) 
which implies ¢,, 0; If oY > 0, then 

i (1) bs 


This argument also suggests that ¢,, * 0. That is to say, the reaction for ¢,, will 
proceed extremely slowly. 

In order to model the ongoing reaction, we obviously have two choices to make 
modifications. One choice is to consider the effect of transport by advection and 


diffusion. The other is to modify the rate laws. 


Effect of compactional flow 


Firstly, let us consider the effect of transport by purely compactional flow. From the 
compaction analysis in the earlier chapters, we understand that u! & 0,0u!/dz = 0 


for the case \ << 1 (slow compaction). Naturally, the effect of compactional flow 


CHAPTER 7. DIAGENESIS: DISSOLUTION AND PRECIPITATION MODEL 116 


is negligible for the reaction. For the cases of 4 = O(1) and X’ >> 1, we have 
u! = O(1), Ou! /Oz = O(1), then O(dxu!)/Oz = O(1). Therefore, Michaelis-Menten’s 
approximation is still valid for ¢x, which means we still have dy & @,. In other 
words we say that the effect of compactional flow is also negligible for the transport 
of species |X]. Smectite moves (at the speed of u*) with the other solids in the matrix, 
and its reaction rate is nearly zero. 

The effect. of transport is only possibly important only if O(¢xu!)/Oz >> 1. This 
implies that the compactional flux should be extremely high. An extremely high 
compactional flux is very rare in natural sedimentation environments and can only 
possibly be generated under very special conditions. In fact, Bjorlykke & Egeberg 
(1993) studied the transport of silica in quartz cementation and concluded that the 
advective transport is not noticeably important in sedimentary basins. 

Therefore, the effect of advective transport is negligible in normal sedimentary 


basins. 


Effect of diffusion 


Since the effect of advective transport is not important, we can simply neglect the 


advection terms in the model equations, and we have 


bm = —Rikidm(ds ar Px) +5 Caras) 


‘ a 7 - 
Obx 1 O°dbx _ 1 [Rikibm(bs — Ox)4 — Rokolox — bs)4)]; (7.78) 


Ot Pe Oz? 


where we have assumed that D = 1. We can easily see that if Pe >> 1, then 
the diffusion is naturally negligible. In the case Pe = O(1), Michaelis-Menten’s 
hypothesis still applies, thus we have dx & @¢,. The effect of transport by diffusion is 
still not important. 

Diffusion will possibly be important when Pe = m,d/Dop << 1. This can be true 
either in a very slow sedimentation environment (small m,) or in fast diffusion process 
(large Do) or in short-distance exchanges (small d). To be more precise, as seen from 


equation (7.78), Pe = O(e) or Pe <<e. 
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In the case of Pe << € (ie. Pe/e << 1), equation (7.78) implies 


Pox 
Oz2 5 


(7.79) 


which means ¢x is a linear function of z. In a pseudo-steady state, the slope is 
determined by mass conservation. If long-distance diffusion dominates the transport 
process, then smectite dissolves in the lower region and illite precipitates in the upper 
region. Therefore, illite should exist even at the top region near the surface, but this 
is contrary to the field observations (illite is rarely found within a depth of 1 ~ 2 
km). This contradiction suggests that only very short-distance diffusion is important 


in natural systems as is suggested by Ahn & Peacor (1986). 


7.5.3 Dissolution controlled or precipitation controlled 


We have seen in the previous discussion that the effect of transport on the reaction 
is negligible. In order to modify the model to mimic the more realistic dissolution- 
precipitation mechanism, we can assume that the solubility of smectite (dissolution) 
is different from that of illite (precipitation) (6 ~ 6/120), and we use the following 


rate functions 
illite 


7, =odml—ox)4, T=(Ox-A4, OF sarc 220 221), (7.80) 
eq 


then the first two equations (7.66) and (7.67) become 


bm = —Rikidm{l — bx}+, (7.81) 
eb = (1 — a) [Rikigm{1 — ox}4 — Roko{ox — 9}4], (7.82) 
Rise 8) Rab See, (7.83) 
with initial conditons 
Om =m and ox = dX. (7.84) 


Since € << 1, the above equations are similar to the model equations for enzyme 
dynamics, then Michaelis-Menten’s pseudo-steady state approximation is valid. This 
fast asymptotics implies Ryki71 & Rokef.. That is 
Rok 


— p(G2—81)O+(F19c,1—829¢,2) 
—e ; ; 7.85 
Rs (7.85) 


Om(1l — dx)4 = A(dx — 6), with A= 
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Solving this equation for dx, we can easily obtain 


(ox — 8), = Se? (7.86) 


Substituting these two expressions into the ¢,,, ; equations, we have 


6-6.,1) PmA(1 — 9) 


5 ott 
Re a 


(7.87) 


J. = (1 = aghem"o-eonrtm 8) — (1 — gg) ferie-eun PAO (799) 


Adding these two equations (7.87) and (7.88), we have 


di +(1—ao)bm =O or oj = (1—ap)(GP, — bm). (7.89) 


Therefore, we only need to solve the first equation (7.87), but it is a nonlinear equation 
whose solution can only be written down implicitly as a quadrature although its 
numerical solution is easily calculated. If A is independent of t, then we can write 
down the solution explicitly. From a geological point of view, we are more interested 


in the following specific cases. 


Equal reaction rates 


If dissolution and precipitation have the same reaction rates with the same activation 
energy (0; = (60, O.1 = 9.2), then A = 1. If the further simplifications 6 << 


1, dm << 1 are used, then we have approximately 
dm = —efO-S cag 


bi = (1— ape OP) b,, (7.90) 


which is the case we discussed before in the dehydration model. We see that the 
dehydration model is a very special case of the dissolution-precipitation model with 
equal reaction rates and very low illite solubility and small volume fraction of smectite 


involved in the transformation process. 
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Dissolution controlled 


When dissolution is the rate limiting process, this is equivalently to two special cases: 
1) By > Bo (Ey > Ey) with Oni y O.2 y 6.2 or 2) O44 > O.2 (Ry < R2) with 
By & Bo & B. 


1) By 2 Bo (9.1 x O.2 x OQ.) 


In this case, A << 1 since 0, = O(1), we then have approximately 
dm(l — ox)4 = AC - 9), 


Ox =H) e= (1-0). (7.91) 
The equations become 


bm = —eO-PiSe (4 = 0). (7.92) 


With a linear temperature approximation 0 = h(t) — z, the solution can be easily 
obtained. 


bm = $9,exp[—-— a? ehalt(t)-21-12e) (7.93) 
Byh(t) on, 


Here we have used the approximation exp(—(,9,) << 1. This solution is obtained 
from an approximation from the original equation (7.55) rather than directly from 
(7.92) since equation (7.92) is only valid in the top part of the region, but the solutions 
we obtained here hold approximately in the entire region. The purpose of writing 
down equation (7.92) is just to show that the reaction rate is nearly independent 
of dm in its region of validity. We can see that dissolution-precipitation will not be 


switched on until a higher critical temperature (O* = 3,0,/(2 > ©-) is achieved. 


2) Oe1 > 0.2 (81 x Bo x B) 


In this case, A >> 1, we have approximately 


bm(l _ bx) + = bm(1 = 0), 


(ox —9)4 = a (7.94) 


The equations become 


bm = —eF(O-Pe) 6. (1 — 8). (7.95) 
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Its solution can be written explicitly if a linear approximation O = h(t) — z for 


temperature is assumed. 


_— 20 A=) B(h-z—®¢,1) 
dm = bm€XP| “Bb e |. (7.96) 


Similarly, the dissolution and precipitation are simultaneously switched on at the 


higher critical temperature 0... 


Precipitation controlled 


When precipitation is the rate limiting process, this is equivalently to 1) 3, < By (Ey, < 
FE) with Oc1 y 0.9 y ie tee or 2) Oe < O.2 (Ry > R2) with By y Bo y ie 
In the former case, A >> 1. With a similar discussion as in the dissolution 


controlled case, we have 
de, = —e®-Pr0014 (4 — 8), (7.97) 
and its solution (for 0 = h(t) — z) is 


1-0 
a= een eens ) ey ae, (7.98) 
Bh(t) on, 
Similarly, the precipitation increases rapidly at a higher critical temperature 0* = 


B20 -/ Pr > On 
In the latter case, A << 1. We similarly have 


bm = —eP(9-822) (1 — 9), (7.99) 
with 
= tbs 0) B(h-z—@¢,2) 
dm = Pm&XP| Bh(t) e ]. (7.100) 


As easily seen, precipitation suddenly increases after temperature is above the higher 


critical temperature QO. 9. 


Equal solubility 


If smectite and illite have the same solubilities, then 6 = 1. We can clearly see that 
all the above discussed processes reach equilibrium very quickly. The transforma- 


tion ceases within a very short time. From a thermodynamic and kinetic point of 


CHAPTER 7. DIAGENESIS: DISSOLUTION AND PRECIPITATION MODEL 121 


view, smectite is equivalent to illite. This is what we have discussed in the previous 
subsection. Fortunately, this impractical case is not of geological interest. 
Effect of transport 


The above results for the dissolution-precipitation model with @ < 1 do not include 


the effect of transport. If we include the terms of transport, we have 


Obm . O(urbm 7 
7 | a = RikiOmtl — Ox }+; (7.101) 
0 O(u' 1, 50 Ue t= 7 ' 
m , ae Peas me) = _ [Rikidm{1 — ox}+ — Roke{ox — O}4]. 


(7.102) 
If Pe = O(1), then the Michaelis-Menten approximation is still valid. From the 
second equation, we can still have equation (7.85). Thus, the transport of [X] does 
not change ¢x noticeably, and the advective term of ¢,, will only change, increasing 
rather than decreasing, ¢,, when the solid matrix moves down but will not change 


the reaction rate. The effect of transport on the reaction rate is, therefore, negligible. 


7.6 Effect of K* and Al* Activities 


Up to now, we have not investigated the interaction of [K*] from K-feldspar. For 
convenience in discussing [At] influence without Al*, the quartz dissolution and 
precipitation processes are not included at the moment. Then the related governing 


equations become (without transport ) 


om = —Rikif, (7.103) 

ox = Riki? — fRokor2, (7.104) 
ob; = fRakga2ho, (7.105) 

oK = = [-asRakiars + asRgk3F3], (7.106) 


by = —agRgks3Fs, (7.107) 
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where € = ox, €x = $x are the equilibrium solubilities. According to Lasaga (1981, 
1984), Ortoleva (1992) and Steefel & Cappellen (1990), we use the following functions 


for the rate laws 


71 = bm(l — $x), (7.108) 
Fo = dnb — 6", (7.109) 
T3 = os(K — Or), (7.110) 


where n is a constant which is usually 1 or 2. « = O(1) is the K-feldspar solubility 
ratio with respect to some reference solubility (e.g. amorphous silica). 
The fast aqueous reaction approximation (€ << 1,€x% << 1) from equations 


(7.104) and (7.106) suggests that 


bm(1— ox) = Abn G% — 6"), (7.111) 
bxGh — 0" = Bo,(k — ox), (7.112) 
with 7 - 
A= = = ae (7.113) 


We see that the above two algebraic equations can be solved for ¢x and ¢x in terms 
of A and B. To understand the [K*] influence more clearly, we are more interested 
in two extreme cases: B >> 1 if K-feldspar dissolves very rapidly while B << 1 if it 
dissolves very slowly. In the following discussions, without loss of generality, we can 
take A = O(1) if smectite dissolution and illite precipitation reactions proceed at the 


nearly same rate. 


7.6.1 K-feldspar dissolution controlled 


In this case, B << 1. By solving ¢x, @x from equations (7.111) and (7.112), the 


governing equations become approximately 
Om = —Rikib; ABR, (7.114) 


0; = fRakza26/Br, (7115) 
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obs = —agR3k36;(1— Bos —0")x, n= 1,2. (7.116) 


We see that the dissolution rate of K-feldspar is nearly a constant (B << 1), but 
smectite dissolution and illite precipitation are very slow as they are controlled by 


the K-feldspar dissolution. 
7.6.2 Fast K-feldspar reaction 
In this case, B >> 1. From equation (7.112), we have 
OK &K. (7.117) 
From equation (7.111), we can easily obtain 
Pm 


= 1—6)—"— for n=1: 


(bm + 20)? + AK m(1 — 8) — bm 
2AK 


Substituting these relations into the governing equations for ¢m, di, Pf, we have (for 


ti 1) 


ox = for n= 2. (7.118) 


: A 
jy = fRakabmaa(l —) 7" — = asRikidm( — Bz, (7.120) 
oy © 0. (7.121) 


By comparing with equations (7.87) and (7.88), we see that the above equations are 
identical to equations (7.87) and (7.88) if we replace A with Ax. Following similar 
procedures, we will have all the results as before. Therefore, the two-step model 
is a very good approximation for the case of fast K-feldspar dissolution and [A+] 
concentration remains nearly at equilibrium solubility. The equations are not much 


different for the case of n = 2. 


7.6.3 Al* activity 


From the model formulations, it is clearly seen that Alt always appears with K*. 


Both cations have similar roles. In other words, they are mathematically equivalent 
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and satisfy similar equations though their chemical roles are different. This means 
that we can always take these two cations as a single combination [K+ — Alt]. There- 
fore, the effect of Alt is essentially the same as K+ and their combination [Kt — Al*]. 


Similar mathematical argument will give similar results as before if we replace x, ¢K 


by bat, bat 


7.7 Quartz Precipitation 


7.7.1 Quartz precipitation controlled 


To include the process of quartz precipitation from the excess silica released by K- 


feldspar dissolution, we still use the same rate laws for 7; and 72, but we use 


P3 = bf(KOs:i — dK osi), (7.122) 
Pa = Osi — Osi, (7.123) 


where 05; < 1 is the ratio of quartz solubility to amorphous silica solubility. Here, 
we have ignored the dissolution process of quartz (7_4 = 0) or can take 74 as the net 
rate of quartz precipitation and dissolution. In addition to the model equations in 
the subsection with KT activity, we still have (without transport) 

‘ 1 >} > a3 a od 

Osi = —|sR3k3r3 = Rakara + Roker], (7.124) 

ES% a4 

where €9; = ¢s; << 1 is the equilibrium solubility for quartz. The Michaelis-Menten 


pseudo-steady steady approximation for dx and @¢s; implies that 


SR3k3r3 + —Raka?s ~ Rakara, (7.125) 
4 
asReksF5 ~ a3RokoFo. (7.126) 
Eliminating 72, we have 
: Raka 
bs (KOs: — OK bs) (ds: — Osi) wi (s + 1)Rakg ( ) 
or 
# CoO iA 
P35 b (KOs; <= OK si) = OF 2 ( Ox) (7.128) 


C+ osoK 
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It is clearly seen that if quartz precipitation is very fast (C >> 1), then og; & 05;, 
the rate law of K-feldspar dissolution is 73 © @f05;(K — @x«), we can then have all the 
similar results as before as in the cases of slow or fast K-feldspar dissolutions. On the 
other hand, if quartz precipitation is very slow (C << 1), then K-feldspar dissolution 
rate is 73 © COs5;(K/¢~% — 1) which is clearly controlled by the rate of quartz precip- 
itation. This is essentially similar to the case of slow K-feldspar dissolution. In the 
top region where the temperature is relatively low, the rate of quartz precipitation 
is very slow. If the temperature increases to some critical value during continuous 
burial, then the rate of quartz precipitation increases dramatically, and this in turn 
switches on K-feldspar dissolution to provide enough A* for illite precipitation, and 
the process of smectite dissolution and illite precipitation will proceed until the reac- 
tion is completed. This reaction series is in line with the recent work by Abercromie, 


Hutcheon, Bloch & Caritat (1994). 


7.7.2 Production of quartz 


The precipitation of extra silica as quartz will have an important effect on porosity 
modifications and reservoir impairments. The calculation of the amount of quartz 
production is obviously needed. For convenience, we neglect the effect of transport. 
By using the full model equations and the pseudo-steady state approximations for 
ox, si,0K, we can easily obtain the relations among 72,73, 74, then we can relate bg 


with bin or di. We have 


‘ s+1 Jot 
bg = -—- bm, (7.129) 
ig 
and 
; stl ‘ 
= aii, 7.130 
Pq fas aie ( ) 
where 


= a7 a3 PmM, 
eae eee ee a 


From the first equation (7.129), we have 


7 ; Cee (7.132) 


dq = 
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Similarly, we can have the expression for amount of K-feldspar consumed 


ae 
ee ne 7.133 
oe fue ee o;Mon ( ) 


or 

gn — —“an( 0, — dn) (7.134) 
With the values of M,, = 367, M, = 60, s = 3, f = 2, Pm/pq © 1,2, = 0.2 (20%), 
then ¢, © 0.065 ( or 6.5%), persed = —0.075 (or — 7.5%) after the completion 


of the S-I reaction. Such a large amount of quartz may have important effect on 


reservoir quality. 


7.8 Summary 


From the above discussions, we see that the reaction-transport dissolution- precipita- 
tion model of diagenesis can reproduce many essential features of the smectite-to-illite 
process if the appropriate reaction rate laws are used based on the known physics 
and chemistry from experimental studies. The detailed investigation of the two-step 
model shows that smectite-to-illite reaction occurs within a narrow region, diagenetic 
window, at a depth nearly ©,, and the reaction processes do not noticeably depend 
on the absolute age of the burial series. More important than geological age is the 
temperature distribution. Long-distance mass transport is negligible in the progress 
of the whole diagenetic process. In addition, we see that the first-order dehydration 
model of diagenesis is a good approximation in the sense of describing the extent of 
progress of the overall smectite-to-illite transformation without much concern for its 
detailed geochemical features. All the results in the case of dehydration model are 
already presented in Chapter 6. 

The full investigations of the whole model in different possible cases reveal that 
K*, Al* cations provided from the dissolution of K-feldspar are very important in 
controlling the progress of the diagenetic reaction. The similarities between the 
present model and Michaelis-Menten’s theory of enzyme kinetics suggest that these 
cations play a role partially like a catalyst during diagenesis. In the case of fast K- 


feldspar dissolution, the two-step model is a very good approximation in describing 
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the whole process. In the case of slow K-feldspar dissolution, the whole process is 
controlled by the rate of K-feldspar dissolution. If there is not enough K-feldspar 
available, then diagenesis can not be completed and may cease in the intermediate 
stage. 

Quartz activity is also a very important factor in controlling the progress of the 
diagenetic reactions. The case of fast quartz precipitation shows no noticeable differ- 
ence from those of slow or fast K-feldspar dissolutions. But slow quartz precipitation 
will hinder the diagenetic process. In the shallow region, the temperature is rela- 
tive low, the rate of quartz precipitation is extremely slow. As burial continues, the 
temperature increases, and the rate of quartz precipitation increases dramatically at 
some critical value of temperature, and this in turn switches on K-feldspar dissolu- 
tion to provide Kt for illite precipitation, thus the process of smectite dissolution 
and illite precipitation will proceed until the reaction is completed. The amount of 
quartz product during diagenesis will possibly have an important effect on reservoir 


quality. 


Chapter 8 


Pressure Solution Creep and 


Viscous Compaction 


The diagenetic modelling in the previous chapter is a transport-reaction model whose 
reaction rate laws do not include the effect of intergranular stress. Obviously, a more 
realistic model should reflect the complexity of the stressed rock system. Pressure 
solution/dissolution has been considered as an important process in deformation and 
porosity change during compaction and diagenesis in sedimentary rocks (Angevine & 
Turcotte, 1983; Tada & Siever 1989). Pressure solution refers to a process by which 
grains dissolve at intergranular contacts under nonhydrostatic stress and reprecipitate 
in pore spaces, thus resulting in compaction. The solubility of minerals increases with 
increasing effective stress at grain contacts. Pressure dissolution at grain contacts is 
therefore a compactional response of the sediment during burial in an attempt to 
increase the grain contact area so as to distribute the effective stress over a larger 
surface. The typical forms of pressure solution are intergranular pressure solution 
(IPS) which occurs at individual grain contacts (Tada & Siever 1989) and free face 
pressure solution (FFPS) which occurs at the face in contact with the pore fluid 
(Ortoleva 1994), but most studies have concentrated on the former one (IPS). In 
spite of its geological importance, the mechanism leading to pressure solution is still 


poorly understood. 
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8.1 Mechanism of Pressure Solution 


Angevine and Turcotte (1983) presented a wonderful theoretical model to study the 
role of pressure solution in the porosity reduction of quartz arenites and the effects of 
grain size, sedimentation rate, and the thermal gradient. Augevine & Turcotte’s work 
was later extended by Birchwood and Turcotte (1994) to give a unified approach to 
geopressuring, low-permeability zone formation, and secondary porosity generation 
in sedimentary basins. A comprehensive review on models of pressure solution was 
given by Tada & Siever (1989). Birchwood & Turcotte (1994) presented more recently 
a brief review on this research subject. 

Extensive studies from petrographic, field and experimental evidence suggest that 
pressure solution is a very complicated process controlled by many factors (Augevine 


& Turcotte, 1983; Tada & Siever 1989; Birchwood & Turcotte, 1994). 


e The main factors controlling pressure solution are temperature, pressure, time, 
grain size and geometry, grain mineralogy, cementation, and solution chemistry. 
These factors do not seem to be simple controlling factors, and may interact 
with each other. The rate generally increases with increase of temperature or 


effective pressure and decrease of grain size. 


e Pressure solution usually involves three successive steps: pressure-enhanced dis- 
solution, diffusive transfer and reprecipitation. The rate of pressure solution is 
thus controlled by the slowest of the three steps, diffusion-controlled or reaction- 
controlled. Most of the existing models of pressure solution assume that diffu- 
sive transport is the rate-limiting step (Weyl, 1959; Coble, 1963; Rutter, 1976; 
Augevine & Turcotte, 1983). 


e The driving force for pressure solution is possibly the gradient of chemical po- 
tential, existing between dissolution and reprecipitation areas, which depends 
mainly on the difference of normal stress, elastic, plastic and surface energies. 
De Boer (1977) has shown that the effect of elastic strains on the gradient of 


chemical potential is negligible. 
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e Pressure solution can start at depths as shallow as 900m in sandstones (Wilson 


& Sibley, 1978) but occurs more commonly at depths of 1 ~ 2 km. 


e Pressure solution is an effective compaction mechanism as well as a source of 
cementing material, especially in sandstones and carbonate packstones. Mass 
transport and redistribution during pressure solution may possibly occur over 


the order of km (Tada & Siever 1989). 


e Pressure solution is probably a combination of plastic deformation and free-face 
pressure solution within and at the edge of the grain contacts (Tada & Siever 
1989). Mass transfer may be carried out by grain boundary diffusion (locally) 
or by bulk diffusion (regionally), depending on the distance of mass transfer. 
Gratier (1984) proposed various pressure solution creep laws for these cases. 
Drewers & Ortoleva (1990) considered pressure solution as a diffusion-reaction 


creep mechanism. 


In spite of all these extensive studies, the operating mechanism of pressure solution 
and the role of plastic versus elastic strain energy as a driving force are still under 
discussion. According to Mullis (1992), two main mechanisms of pressure solution 
are possible. One mechanism assumes that increased solubility at the grain boundary 
sets up a concentration gradient resulting in mass transfer by diffusion into the pore 
spaces. The diffusive transfer could occur by thin water film diffusion adsorbed to the 
grain boundary (Weyl, 1959; Rutter, 1976; Augevine & Turcotte, 1983), or by bulk 
diffusion or through fluid ‘island channels’ (Raj & Chyung, 1981). An alternative 
mechanism is the undercutting model which supposes that the increased solubility 
at the grain contacts results in preferential dissolution at the rim of grain contacts 
leading to undercutting and brittle failure (Bathurst, 1958) or plastic deformation 
(Pharr & Ashby, 1983; Tada, Maliva & Siever, 1987; Pytte & Reynolds, 1989). 

Many experiments have been carried out to investigate the mechanism of pressure 
solution, but no evidence has been found for grain undercutting though neither is 
direct evidence for adsorbed thin water films convincing. It is worth pointing out that 


the theory of (adsorbed thin) water film diffusion (WFD) is theoretically favoured by 
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the concept of very high disjoining pressure (~ 270 MPa) of removing a monolayer 
from a thin film (De Boer, 1977; Rutter, 1976; Tada, Maliva & Siever, 1987; Mullis, 
1992). However, some controversy still exists about which of the two mechanisms 
is more appropriate, as they produce the same creep rate in the simplest case of 
diffusion-controlled creep; the only difference lies in the interpretation and values of 


parameters such as the effective grain-boundary diffusion coefficient. 


8.2 Mathematical Model 


For the convenience of investigating the effect of pure pressure solution, we will begin 
by neglecting diagenetic reactions such as the smectite-to-illite transformation and 
will assume a single species only such as quartz. Extensions will then be made to 
reactive multiple species during diagenesis. From the discussion in the previous sub- 
section, we will also assume that the dissolution-diffusion-precipitation process only 
occurs locally on a grain scale. Based on the existing models of pressure solution 
(Weyl, 1959; Rutter, 1976; Angevine & Turcotte, 1983; Nielsen, 1986; Mullis, 1991) 
and models of compaction (Fowler, 1990; Stevenson & Scott 1991; Audet & Fowler, 
1992; Birchwood & Turcotte, 1994), the present model is written as 


Conservation of mass 


(1-4) +¥-[(1—¢)u'] =0, (8.1 
2 9d) =0, 2) 
Darcy’s law 
ou! —u’) = -F (Vp! + ai) (8.3) 
Force balance 
V-o* —V[(1—a)p'] — p9j = 9, (8.4) 


where € is bulk viscosity, a. is the effective stress, p, is the effective pressure, j is the 
unit vector pointing vertically upwards, k is the matrix permeability, py is the liquid 


viscosity and p! is the pore pressure. 
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A rheological relation and a constitutive law are needed to complete this model. 
If an Athy-type law is used to relate effective pressure p, with porosity @, we will 
again return to the model we have analyzed in the earlier chapters. To emphasise the 


features of pressure solution creep, a different type of constitutive law is expected. 


8.2.1 Constitutive creep laws 


There are generally two ways to make mathematical formulations of creep laws for 
pressure solution. One way is to derive creep rate in terms of concentrations, grain 
size and geometry (usually spherical or cylindrical packings), effective stress, grain 
boundary diffusion. Most models fall into this category (Weyl, 1959; Paterson, 1973; 
Rutter, 1976; Angevine & Turcotte, 1983; Mullis, 1991, 1992; Shimizu, 1995; Lehner, 
1995; Schneider et al, 1996). This allows us to include the detailed reaction-transport 
process in a simplified relation between strain rate and effective stress although fur- 
ther simplifications are usually assumed such as steady-state dissolution and local 
reprecipitation along the grain boundary. An alternative method of formulation is 
simply to assume a viscous law, as is done, for example, in modelling magma transport 
(McKenzie, 1984; Fowler, 1990). The latter treatment does not describe the details 
of the pore scale dynamics. The connection between these two kinds of formulations 
can be easily seen from the constitutive law used in the formulations. 

In the formulations of the first kind, the Weyl-Rutter creep law is widely used 
(Weyl, 1959; Rutter, 1976; Angevine & Turcotte, 1983) 


; Arco wD 
C= ys (8.5) 


where o is the effective normal stress across the grain contacts, A, is a constant, Co 
is the equilibrium concentration (of quartz) in pore fluid, p, d are the density and 
(averaged) grain diameter (of quartz). Dyy is the diffusivity of the solute in water 
along grain boundaries with a thickness w. This relation holds for the case of steady- 
state diffusion and no grain-boundary sliding. 

The relation between porosity ¢ and volume strain e depends the grain geometry 


and packing texture. Weyl (1959) and Augevine & Turcotte (1983) used the following 
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relation 
1 
615 (8, 8.6 
for regularly 3-D spherical packing, while Schutjens (1991) used 
_ go-o 
Oe 5° (8.7) 


in explaining the experimental compaction of quartz. D,y also varies with temperature 


T. In fact, Augevine & Turcotte (1983) wrote D,,(T) as 
Fea 
Dy(T) = Dye FT , (8.8) 


where Hq is the effective activation energy with a value of 3 ~ 6 kJ/mole or even 
much lower (Angevine & Turcotte, 1983; Nakashima, 1995; Shimizu, 1995). From 
the values of the diffusion coefficient used by Paterson (1995) in quartz-water and 
rocksalt-water systems at 300, 600, 1200 K, we get an estimate value of H.g = 0.65 
kcal/mole. Rutter (1976) and Angevine & Turcotte (1983) pointed out that these 
values are only estimations. 

In the classical formulations, the following constitutive laws are often used (Roberts 


& Tabor, 1971; Paterson, 1973; Mullis, 1991) 


o o 
RY? and cage (8.9) 


Cc = coexp(— 
where Wo, 09 are constants depending on the properties of the thin film, and v is the 
molar volume (of quartz). These constitutive laws, though experimentally based, are 
essentially theoretical simplifications as in the case of Athy’s law p. = p.(@). 


In addition, the diffusion coefficient D,p also depends on the porosity ¢. According 
to Archie (1942) and Paterson (1995), we have the following Archie’s law 


Dov = (2)"D,, (8.10) 


o 
where Dp is the diffusion coefficient at the initial porosity @9, and n is the exponent in 
Archie’s law. The value of n has been determined empirically to be 1.3 for uncemented 
sand-like granular media and 2 for a wide range of cemented rocks (Paterson, 1995). 
Experimental studies show that A; in the creep law (8.5) is not a constant and 


depends on temperature T (Raj & Chyung, 1981; Augevine & Turcotte, 1983; Spiers 


CHAPTER 8. PRESSURE SOLUTION CREEP AND VISCOUS COMPACTION 134 


& Schutjens, 1990; Dewers & Hajash, 1995). To incoorporate this, (8.5) is usually 


modified as 
Aalo wDgp 
RT 


where R is the gas constant, and Ay = av?,/vy,0. a is a factor depending on the 


é€=A(¢,T)o, with A($,T) = (8.11) 


grain geometry and stress distribution. 1, is the molar volume (of quartz) and vy,0 
is the molar volume of water. The constitutive relation (8.11) is originally for the 
one-dimensional case. We can extend it to a more general form by writing 

RT@ 
—————— €xk. 8.12 
Agco wD, oe 
Note that p. = —o and &, = V-u®. With this, (8.12) becomes 


RTE 
fp pa. 8.13 
ae eerie a 


CO = 


which is equivalent to the following compaction law 
De = —EV.U*. (8.14) 


This was first used by Birchwood and Turcotte (1994) to study pressure solution in 
sedimentary basins by presenting a unified approach to geopressuring, low perme- 
ability zone formation and secondary porosity generation. The compaction law is 
analogous to dislocation creep controlled viscous compaction laws used in studies of 
magma transport in the Earth’s mantle (McKenzie 1984, Fowler 1990). 

Another way of formulating constitutive laws for pressure solution creep is to con- 
sider it as a viscous compaction mechanism. The creeping process under effective 


pressure p, can be formulated as 


d(1— ¢) 
dt, 


where d/dt is the material derivative 0/Ot+ uS-V, following the solid matrix. Rewrit- 


= K($,T)pe, (8.15) 


ing the equation (8.1) of mass conservation 


= aie —(1-¢)V.u%, (8.16) 
then (8.15) becomes 
1-@ 
Pe = — V.u’, (8.17) 
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which is equivalent to (8.13). Therefore, we can say that the two kinds of formulations 


can be unified as a single creep law such as (8.14). 


8.2.2 Derivation of creep law 


The approach of deriving the law of pressure solution creep depends on the underlying 
mechanism. The classical theoretical consideration used by Weyl (1959) and Rutter 
(1976) assumed a grain-boundary diffusion film of constant thickness and diffusivity, 
while others used the concept of a roughened, fluid-invaded non-equilibrium contact 
structure (Raj, 1982; Lehner, 1990; Spier & Schutjens, 1990; Lehner, 1995). Shimizu 
(1995) presented a kinetic approach extending Coble’s classical treatment of grain 
boundary diffusion creep (Coble, 1963) by including the kinetics of quartz dissolu- 
tion/ precipitation reaction. Shimuzu’s (1995) derivation is instructive although the 
boundary conditions used in his 1-D diffusion model are questionable. This 1-D ap- 
proximation is only valid for a closed system when the thickness w of the water film 


is small with respect to the grain diameter (d). 


grain matrix 


ee 


—>S  <— 
<< 


a 
| 


Figure 8.1 Water film diffusion model of pressure solution creep in which 
dissolution occurs in the contact region and reprecipitation takes place 


along grain boundaries in pore space. o* is the effective normal stress. 


Existing pressure solution creep laws can produce some typical features such as 


pressure and grain-size dependence of creep rate, but have essentially remained re- 


CHAPTER 8. PRESSURE SOLUTION CREEP AND VISCOUS COMPACTION 136 


stricted to macroscopically closed systems, with negligible long-distance transport 
in the pore fluid, and they are somewhat biased toward grain-boundary diffusion- 
controlled pressure solution creep. This shortcoming of the creep laws eliminates the 
effect of solution transfer over large distances. The consequence is that the coupling 
effects between solute transport and pressure solution deformation under possible 
open system conditions encountered during sediment diagenesis or metamorphism 
have received limited attention (Dewers & Ortoleva, 1990; Lehner, 1995). Augevine 
& Turcotte (1983) pointed out that future quantitative modelling of sediment dia- 
genesis should incorporate temporal variations in subsidence rates, spatial variations 
in lithology, and heat flow (Turcotte & Schubert, 1982). Lehner (1995) investigated, 
for the first time, the creep law of pressure solution in open fluid-rock systems by 
recasting the classical Weyl-Rutter model of intergranular pressure solution in terms 
of a (linear) phenomenological creep rate law. However, Lehner (1995) left an open 
question concerning the validity of the postulated simple creep law due to the uncer- 
tainty of the assumed phenomenological rate constant K%, and suggested that a new 
generation of diagenetic models should describe an open system including the effect 
of composition, fluid pressure, temperature as well as solid stress state. 

Now let us consider the intergranular contact region as a disk with a radius r = L. 
Let J(r) be the radial component of solute mass flux, é be the average strain rate, 
and v is the uniform shortening velocity of the upper grain relative to the lower grain 
due to the pressure solution creep. The kinetic relation between v and é€ is (Lehner, 
1995) 

v = éd. (8.18) 


For simplicity, we assume that the film thickness w is constant and the diffusion is 


near steady-state as Rutter (1976) and Lehner (1995) did. Mass conservation gives 
QrrJ(r) + psrr?v = 0, (8.19) 
where the flux J(r) obeys Fick’s Law 


J(r) = —Daw (8.20) 


dr’ 
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The steady-state solution of concentration c(r) for the boundary condition c, = 0 at 


r=0,c=caatr=Lis 


ofr) =e — 7 Dae — pr), (8.21) 
g 


The parabolic change of concentration c(r) implies that the stress o(r) should be 
heterogeneously distributed in the contact region. From the relation (8.9), we have 


afte) 


o(r) = as (8.22) 


where we have used the condition o°(r) = 0 at r = L. Let o be the averaged effective 


stress, then 


L 
tL?o = i, 2ro°(r)rdr. (8.23) 
0 
Combining (8.22) and (8.23), we have 


2QRT tl psed 
Se) InsP a dirs 8.24 
7 mi el reese Pale cep 


Using (8.18) and integrating by parts, we have 


2 
= ——_|(1 —- —~ — — 2 
oO ie [(1 Bp ad BL*) — 11, (8.25) 
where 7 
psed 
— ; 8.26 
AcoD yw ( ) 


By defining a critical effective stress a, (and equivalently a critical creep rate é, ) 


when BL? = 1 


RT, AcoDgyw 
e Vm i psL?d ( ) 
(8.25) can be rewritten as 
or € € 
—=f1-(1- =)Infi - —)]. 8.28 
2 = (1-1-2) 5) (8.28) 


A typical value of o, is about 95 MPa with values of T ~ 300 K, R ~ 8.31 J mol"! 
K-t, and vn, ~ 2.6 x 107° m3 mol™?. 
Clearly, if |7|< o-, we have 


~ AUmcoDgw _ _ 16UmcoD gow 


= ?gb f 8.29 
RTpdi?. ~~ = RT og ony) 
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which is exactly the creep law (8.11). Here we have used L = d/2. A different choice 


of L = O(d) will only introduce an additional shape factor into the above relation. 
Under upper-crustal stress conditions ¢ < 100 MPa (Zoback et al, 1993), the above 
approximation is valid as we expected. At higher stress states, we can use |a|>> a, 


then (8.28) becomes 
i Aco Dgpw 
pad L? 
Let L? = 4d?/a,, and a, = O(1) is a shape factor. The above relation (8.30) becomes 


[l —e7 Fr]. (8.30) 


a AsCoD gow 
psd® 


Umno 


[l —e7 BF], (8.31) 


which is consistent with Dewers and Hajash’s empirical law derived from a quartz 
compaction experiment (Dewers & Hajash, 1995; Siese & Spiers, 1997). It is worth 
pointing out that the creep law (8.31) degenerates into (8.29) when v,,0/RT < 1, 


but it may be inaccurate when |a|~ o¢. 


8.2.3. Equation of motion 


For an open system, we expect that a source term will be introduced into the macro- 
scopic equation of mass conservation. Now the porous medium consists of two phases, 
the solid matrix (quartz) and the pore fluid (dissolved silica and water). Let ¢,, be 
the volume fraction of the solid, cs; and c,, be molar concentrations of the dissolved 
solid species and water in the pore fluid, respectively, r,, be the rate of mass dissolved 
by pressure solution, and r, and r_ be the rates of dissolution and precipitation (of 
quartz) on free surfaces of grains where the effective pressure 0° = 0. The equations 
of mass conservation now become 


Conservation of mass 


Oem Ca ee Mm 
OE + V.[onu | — ( lm ry F r_) Pos ’ (8.32) 
a — V.(DeVeg: — u'des;) =tm+r+— 1, 338) 


dmto=1. (8.35) 
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In the above equations, we have used volume fractions (¢,, and @) together with molar 
concentrations (cs; and c,). For convenience in the following non-dimensionlization 
and analysis, we will rewrite the above equations solely in terms of volume fractions 
os; (of the dissolved solid species) and ¢,, (of water) to replace the molar concentra- 
tions (cs; and c,). Clearly, 6 = 5; + dy is the porosity and ¢,, = 1— ¢. Thus, we 
have 


Conservation of mass 


Oil - Min 

SOF + VIG = 80!) = (tm 4 +), (8.36) 
oto esi) + V.[(¢ — ¢s:)u'] = 0, (8.37) 

es — V.(DV és; — u'ds:i) = (tm t 74 — 1) . : (8.38) 


where M,,, and Mg; are the molar weights of the solid and the disolved solid species, 
respectively. p,;, and p; are the densities of the solid and the pore fluid, respectively. 
D is the diffusion coefficient in the pore fluid. If the process is only pressure-enhanced 
and there is no free surface reaction involved, then we have only one source term rp. 


The source term 7, is 
Tm = Np,unL? = Np,édL’, (8.39) 


where N is the number of grains per unit volume. Substituting (8.29), we have 


_ Atm CoD gw 


No. A 
Wee RT o (8.40) 


For 3-D packed spherical grains, N = Ag/d?, we finally obtain 


— 4AgT Un coD gow 


Ym = (8.41) 


where Ap = O(1) is a shape factor which is 1 for cube-shaped grain packing and 6/7 
for spherical grain packing. 
From the kinetic theory of the quartz-water system (Rimstidt & Barnes, 1980; 


Paterson, 1995), quartz pressure solution is described approximately by the reaction 


SiO2(solid) + 2H,O(aq) = H4SiOu(aq). (8.42) 
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The free surface dissolution rate r, and precipitation rate r_ depend on the temper- 


ature, pore pressure p and the concentrations. 
ra = kGlL, p) Osis. G6; (8.43) 


t= KhALD) Grids, (8.44) 


where a is activity. For free surface dissolution and precipitation in dilute solutions, 
we may assume 

ano © 1, asio, 1, Anysio, © Psi- (8.45) 
Under the natural conditions of sedimentation, the dissolution rate constant k+ 
and precipitation rate constant k_ change with temperature 7. Based on Rim- 
stidt and Barnes’ (1980) theory, we have ky ~ 2.0 x 10-°molm?s~! (at 300 K) 
to 59 x 10-® molm’s~! (at 600 K); k_ ~ 1.5 x 10-4 molm?s~! (at 300 K) to 8.0 x 
10-4 mol m?s~! (at 600 K). 


8.2.4 Compaction relation 


Now the total strain rate €;; 


1 dug  Ous 
= oA ; 8.46 
5 7 2 ( Ox; Ox; ( ) 
in the sediments is considered to be partly elastic ¢{; and partly viscous €?; 
big = Ef + Es (8.47) 
whose Maxwell formulations are in the following form 
Vu . Be lastic) — y(@, T)pe(viscous) (8.48) 
w= elastic) — 7(¢,T)pe : 
(1 — )p.(¢) dts 
or 
doy, 
Ekk = g(¢) a — (9, T) Che, (8.49) 


where g(@) is a known function of ¢. We also have used equations (8.16), (8.14) and 


Athy-type law Pe = Pe(Q), p.(¢) <0. 


@ ., WmCoDow T —To. (T-T) 24 
T)= ff ——(]1 — a0 8.50 
(8,7) = (Sy emcee — FT (8.50) 
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where Jo is the surface temperature at the top of basin. Ag is an effective activation 
energy (Augevine & Turcotte, 1983). For convenience in the following discussion, we 
will prescribe the temperature distribution. 

In principle, we can write a generalized Jaumann-type relation in a corotational 
frame of reference for (8.49) by considering the coordinate objectivity (Fowler & Noon 
1995; Khan & Huang 1995). For simplicity and clarity in the following analysis of 
the effect of pressure solution creep, now we will mainly focus on the purely viscous 


compaction and use the following compaction relation 
De = —EV.U"," (8.51) 


which was first used by Birchwood & Turcotte (1994) to present a unified approach 
to geopressuring, low permeability zone formation and secondary porosity generation 


due to pressure solution in sedimentary basins. 


8.3 1-D model and Non-dimensionalization 


8.3.1 1-D model 


For simplicity, we let a = 0 and o = o§ be the averaged effective stress (z—component). 


The 1-D model equations then become 


a-¢)  d-dw] Mm 
at T Be — ( Tm ry Wr re) De ’ (8.52) 
Ag — si) , AlO—=osi)ul] _ 
oy ae = 0, (8.53) 
Oo, OO, Obs: - Msi 
Ot — Oz (D Oz =~ u'dbsi) = eas of Ty 7? i) pI ) (8.54) 
(ul — ut) = "(PP + pig), (8.55) 
Oz 
0. 
5, 7 Pa(b— 8) + ida = 0. (8.56) 
O Ss 
a = (eT ho. (8.57) 


where 7(¢, 7) is a function of porosity ¢ and temperature T. 
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Now we have 6 equations for 6 unknown variables: two for porosity ¢ and ¢5;, 
two for velocities u°, u!, two for effective stress o and pore water pressure p. The 


boundary conditions are 


oO 0, Pp 0, Q 0; Osi = bsi, at z= A(t), (8.58) 
uw=—0, ub =0 at z=), (8.59) 
A(t) =m, +’. (8.60) 


8.3.2 Non-dimensionalization 
We scale the effective pressure —o with (p; — pr)gd 
—o = (ps — pi)gdp, (8.61) 


so that p = O(1)). We will define the length scale d by equation (8.72). We also scale 
pore pressure p with (p, — px)gd, permeability k with ko, time t with d/m,, z with d, 
k4 with k9, k_ with k°, dg; with és; and putting 


T=To+ 0, (8.62) 
then the dimensionless model becomes 
Conservation of mass 
o(1 — O|(1 — d)us sia a : 
( ey 9) ( 9) = —AI7p — a\(Rif, —R_f_), (8.63) 
= : = Yat 
HO — cbs) , A6— etal _ 9 as 
Oos; 1 O ~0bdg; 1 ae ‘3 e 
SS (DSS —ulds)) = ALG + (Rave — RF, (8.65) 
Darcy’s law 
I s T Op 
e(w — us) = Ala +r), (8.66) 
Force balance 
Oy Y 
ao Pe) SU) 0, (8.67) 
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compaction relation 


Ou® 
— = -l4p 8.68 
De AD, (8.68) 
where 
ko(ps — pr) WmCoDow(ps — pr)gd? 
an ao Ps PUT. ya Om CoV OW Ps Pug (8.69) 
pms : RTop,d?ms5 
= (D1 e)e®, F=(Ly", D=e®, 
Po Go 
Faqod god Ar Ag Mnpi 
= A= = . 
Ben RKoTe’ fen Kyi ’ a ’ ai PmM 5; (8 70) 
Igy hts kod R kod pes math __ Al 
Pe mis ’ = is ’ oa ’ pe pL 
3 “ E.qod E_qod 
— °f+9(4 ry nee, eee a alo 71 
ry € ( ¢), = € Osi; Bs KoRT2’ Be KoRT2’ (8 7 ) 


and € = ¢g; ~ 6 x 10-® < 1 is the quartz solubility in the pore fluid. Now we can 


define a length scale d by setting [ = 1 


Ti Ss B ? s 
= FE DSE pe», (8.72) 
AVnCoDow(Ps — pi)g 
The related boundary conditions become 
p = 0, Q = o; OF oo Peis at z= h(t), (8.73) 
wa=u =0;. at-z=0, (8.74) 
A(t)=m.+u’. (8.75) 


8.3.3 Values of parameters 


By using the typical values of p; ~ 10?kgm7?, ps ~ 2.5x10?kgm7?, ky ~ 107!8m?, pp ~ 
10-3>Nsm?, d ~ 10-*m, R ~ 8.31Jmol 1 K7!, v, ~ 2x 107>m?mol ja ~ 
16, Hig ~ 3kcalmol™', To ~ 300K, co ~ 1074*+M, woDgs, ~ 1 x 107% ms? (Rut- 
ter, 1976; Gratz, 1991; Birchwood & Turcotte, 1994; H. Ockendon & J. R. Ock- 
endon, 1995), Go ~ 1x 10°’ Pa,v ~ 0.2,d ~ 900mm, m, ~ 107 ms™!, ko ~ 
10-1 mol m?s71, k° ~ 5x 107! molm?s!, Ey ~ 51.4kJ mol, E_ ~ 34.3kJ mol, 
then we have 


hel, PA, Pe © 30, r= 0.6, (8.76) 
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Bm = 0.2, 3, = 0.1, m = 8, = 2, a, & 0.4, (8.77) 


R. 0.1, Ri AS, Bx 2 1.8, Gow, Aw 0.8, (8.78) 


where we have used M,,, = Mg; and |qo/Ko|= 30° C/km (thermal gradient). 


8.3.4 Effect of transport 
Because ¢ < 1, the pseudo-steady state theory applies, this means that 
AIYp + ay (Rats = R_f_) _ O(e) ~ 0, (8.79) 


which is an algebraic equation for ¢s; in terms of p, 0 and @, i.e. G5; = bs: (D, 9, ¢). 
Adding the first three equations (of mass conservation), we have 


due + o(ul—us)]  € 0 -Obs, 1 
Oz = Peds! Oz )4 o 


~1)0(6), (8.80) 


which implies that the effect of diffusion is only significant in a characteristic diffusion 


€ 
dp = \ 5,4: (8.81) 


which is approximately 0.4m. The diffusion length dp «< d shows that reprecipitation 


length dp 


essentially occurs locally, and the effect of long-distance transport is negligible in the 
natural sedimentation environment. Furthermore, « < 1 implies that ¢ — eds; & @. 
The locality of reprecipitation of the dissolved species by pressure solution enables us 


to model the process by using the reduced equations 


A(1— 4) | Ad — ¢)u*| 


=) 8.82 
Ot Oz ; (eee) 
Og , Agu’) 
— + ==) 8.83 
Ot Oz ; eed) 
Adding these two equations and integrating, we have u’ = —¢(u'—u®) in a barycentric 


frame. Solving pore fluid pressure p from the force balance equation, and substituting 
into Darcy’s law, then we have only three coupled governing equations left in the 


model. They are 
dl — ¢) 
Ot 


Ol(1 — d)u] 
Oz 


= 0, (8.84) 
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s_ _yp°P 4 
ue = ARLE + (1-¢)]. (8.85) 
0 s 
5, = Dib, (8.86) 


which are essentially the equations given by Birchwood & Turcotte (1994). The 


related boundary conditions become 


p=0, d6=¢o, at z= hit), (8.87) 
i= 0y are, (8.88) 

; ~_OD 

A(t) =m, — MS (=e) (8.89) 


8.4 Viscous Compaction 


If we put ==1/T, éf= 1/7, then the governing equations become 


O1— 4) | All—¢)u*] _ 
oy | De =: (8.90) 
ts Op 
°=—)\k[— + (1-4). OL 
w= MSE + (1 9) (8.91) 
= OU 
p= ath aS (8.92) 
where € is porosity-dependent and f is temperature-dependent. i.e. 
is b._ 
a(t 8.93 
c= be? (8.93) 
f=(1+B,0)eFr®, (8.94) 
Eliminating p, the governing equations become 
0g _ AU — )u*) 
Ot Oz ; eee) 
>, O ,~zOu® 
2 NES —(1-4)]. : 
uw = MET (EF) - (1 - 0) (8.96) 
Now the boundary conditions are 
Ou® 
ae 0, d=, at z=A(t), (8.97) 
(iSO. Bh eb, (8.98) 
A(t) =m,+u'. (8.99) 


This problem is very difficult to analyse. In the rest of this chapter, we will mainly 


solve it numerically and give some asymptotic analysis when it is possible to do so. 
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8.5 Numerical Results and Analysis of Viscous Compaction 


From the governing equations and their boundary conditions, we understand that the 
physical model suggests that A, =, k, f are all positive, 6 > 0, > 0 and u® < 0. 


The non-negativeness of p implies 


Ou® 
< il 
ap <0 (8.100) 


which means u*® monotonically decreases as z increases. Thus, u* reaches its minimum 
i= (2 = h) at z= and its masamum. uv? = 0 ate — 0. 

It is seen that the coupling of the temperature O with ¢, u® in the governing 
equations appears in the form (1 + 3,0)exp(—G,,0) which complicates the analysis. 
In order to separate the effect of temperature change from that of viscous compaction 
with a constant temperature distribution, and to compare with the existing results, 
we will mainly discuss viscous compaction without temperature change (3,, = 0, = 0 
or f = 1). These simplifications are in fact reasonable since 6,, ~ 0.2 < 1 and 


8.5.1 Slow compaction \ << 1 with = = O(1) 


The numerical results for the case of small is shown in Fig. 8.2. A boundary layer 
clearly occurs at the basement. 

For the case of \ << 1, we can expect that u®’ = O(A) << 1 since k = O(1), 
this implies that 6 = O(A) or @ & ¢o and thus k = 1, € & 1 which are in line 


with the numerical results. With these simplifications, the model equations become 


approximately 
dt = (1 — do)ué, (8.101) 
AZu;, — uu? = (1 — op). (8.102) 
The outer solutions are ¢ = ¢o, u® = —X(1 — go), and there is a boundary layer at 


the base, for which the effective boundary conditions can be written as 


270, ¢—+¢0, as 2-700, (8.103) 
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The solution for equation (8.102) can be easily written as 
u’ = A(1— do)[e P= — I]. 


Substituting this solution into (8.101), the solution for ¢ is approximately 


b= $y — (1 — $o)4f Ate FE, 


ay 


and 


h=1—AX(1— do). 
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(8.104) 


(8.105) 


(8.106) 


(8.107) 


We see that a boundary layer is developed at the base with a thickness of the order 


of VAS. The comparison with numerical results shows good agreement (Fig. 8.2). 


Dashed: Solution 
o.2b Solid: Numeric 
0.15 
= 
> 
oO 
= 
mo} 
x) 
[so] 
GH 0.1 
0.05; 
5_ 
0 f ———— 
0.3 0.35 0.4 


Porosity 


Figure 8.2 Comparison of solutions (as dashed curves) in the boundary 


layer with numerical results (as solid curves) in the case of slow creep 


(A <1) for t = 2,3,5. Z is the scaled height z/h(t). 
8.5.2 Fast compaction \ >> 1 with = = O(1) 


Numerical Solutions 


The numerical solutions for A = 100 at different times (t = 1,2,3,5,8) are shown in 


Figures 8.3-5. 
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Porosity vs scaled height [Fast Creep] 
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Figure 8.3 Porosity profile versus scaled height Z = z/h(t) at different 
times (¢ = 1, 2,3,5, 8). 


Porosity vs depth [Fast Creep] 
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Figure 8.4 Porosity profile versus depth z — h(t) at different times (¢ = 
12 325;8) 
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From Fig. 8.3 and Fig. 8.4, we can clearly see that the porosity profile is nearly 


in a parabolic shape in the region near the top, and moves as time t increases, which 


suggests that there exists a travelling wave solution in the top region (Fig. 8.4). 


On the other hand, the solutions at longer times suggest a different feature below 


the transition region where the compaction becomes permeability-controlled as the 


porosity decreases so that Ak << 1. 
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Velocity vs height 
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Figure 8.5 Velocity profile versus height h(t) at different times (t = 
1,2, 3, 5,8). 
Figure 8.5 shows the velocity profile at different times, and this profile suggests 
that the velocity at the top tends to be a constant as t increases. In other words, this 


means that h = ms + u® becomes nearly constant. 


Poro-elastic compaction versus viscous creep 


To study the behaviour of creeping compaction, it is helpful to compare numerical 
solutions with the counterparts for poro-elastic compaction. Figure 8.6 shows such a 
comparison with values of \ = 100 and t = 5. 

It is clearly seen that poro-elastic compaction behaves differently from viscous 
compaction in the top region in that the former decreases more rapidly than the 
latter, but they behave in a similar way in the lower region. This is because the high 
permeability near the top will enable poro-elastic compaction to proceed rapidly, 
leading to a nearly exponentially decreasing porosity profile in the top region, but 
the low effective pressure near the surface will only make viscous creep proceed slowly, 
resulting in a nearly parabolic profile of porosity evolution at the top region. On the 
other hand, the porosity decreases to values lower than the critical value ¢* below the 
transition region, where the permeability is low enough (ie. Ak < 1) to retard the 


compaction or creep process, so that both processes essentially become permeability- 
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controlled, resulting in a similar profile in the lower region. 


Poro-elastic and viscous compaction 


Dashed: Poro-elastic i 


0.8+ Solid: Viscous creep A 


Scaled height: Z 
Oo 
oa 
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Figure 8.6 Comparison of poro-elastic compaction (dashed) with vis- 
cous compaction (solid) for \ = 100) at t = 5. Here we have choosed 
the parameter values in equation (8.72) in such a way that the length 
scale defined by (8.72) is equal to the length scale defined by (2.48) in 
the poroelastic compaction. We then can make the comparison with 
the same length scale. The choice of other combination of the parame- 
ter values in (8.72) will make the two curves look very different in this 


figure. 
8.5.3 Analysis for \ >> 1 


From the governing equations (8.95) and (8.96), we see that Ak appears always as a 
combination. Because k decreases dramatically as ¢ decreases, we can expect that 


the value of ¢ when Ak = 1, or 


Ind 


b=b =He™, (8.108) 


will define a transition as in chapter 4. 


Since \k >> 1 holds in the top region, then we have 
1-@x =—(E—-). 8.109 


Substituting this into equation (8.90) and interchanging ¢t and z in differentiations, 
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we have 
0.0, ~0us 0 ~Ou* 
alae’ 


5) tua (ES -)] = 0. (8.110) 


Rewriting the boundary condition uf = 0 at the top by multiplying by € and differ- 


entiating with respect to t, we have 


0 20u° 2 CO .0U? 
ai Oz )4 hae Oz ) 


=0. (8.111) 


Using this relation and integrating equation (8.110), we obtain 


0 ,2du*, 0 ,20u®, __tins(1 — do) 
aon on ee ee 


(8.112) 


fl 


Let p = p, u = u’, and for simplicity choose € = 1, then we have the following 


equations 
pet up. = ms(1 — 0), (8.113) 
p=—Suz. (8.114) 
The boundary conditions are 
p=0 at z=A(t), (8.115) 
and 
C=O" ah 20: (8.116) 


The characteristics of (8.113) are 
z=u and p=™m™,(1— ¢o), (8.117) 


The boundary condition (8.115) can be written as 


(=a 2h ype 0: (8.118) 
By integration, we have 
p=m,(1 — bo)(t — 7), (8.119) 
t 
a | u(s,T)ds + h(r), (8.120) 


and 


ty(T) 
h= -| u(s,7)ds, (8.121) 
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where t;(7) is the time t at the basement z = 0. 
Changing variables from (t, z) to (t’,7) (in fact, t/ = t), we get 
1 t 
is = ie and: 2S is +f u-(s,T)ds. (8.122) 
Zr T 
Now we rewrite equation (8.114) as 


el rz. o) 


—U, 


S eee ts + fru,(s,7)ds’ ete) 
whose integration leads to 
mr a ee oe ms Feucls of} (8.124) 
or 
[ u,(s,T)ds = tnglew Et — lj, (8.125) 


Differentiating with respect to t, then integrating u, with respect to 7, and noticing 


that f,(t — 7) = —f,(t — 7), we obtain 


u=h—me Ee, (8.126) 
Using the boundary condition u = 0 at z = 0, we get an expression for h: 
0=m,[1 — en ED (t5(7)-7)") +(h — ths), (8.127) 
that is 

f= rage EO), (8.128) 

Substituting the solution (8.126) into equation (8.120), we have 

t ; _ ts—do) (g_7)2 : . 

= i {ring[L — en = | + (A(t) — rh) }dt = A(T) — z, (8.129) 


Changing variable t = tT + vera the above equation becomes 


tis (1-99) (t—r) 
ee ae e~”' dv — [h(t) — A(r)] = h(r) - 2. (8.130) 
eles ms( = do) 
a ea) (6.131) 


Here, we see that h — z > eat 5p) a8 t— T — 00. 


—0 
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It is clearly seen that u® and p are only functions of h(t) — z, which implies that 


@ is also a function of h — z. The equation of conservation of mass then becomes 
ho! + (1 — dul’ = 0. (8.132) 


Integrating this equation (8.132) and using the boundary conditions at z = h, we 
have 


hd + (1 — $)u® = hoo + (1 — do)(h — ths), (8.133) 


(1 — )(h—u) = m.(1 — 0), (8.134) 


Substituting solution (8.126) into this equation, we have 
peas Gee a (8.135) 


To find the solution for h(t), we rewrite equation (8.128) in terms of s, = t, —T as 


: : ths(1—¢0) .2 
h=m,e =. 


(8.136) 


When z = 0, equation (8.131) becomes 


Nae Sa ll Mel T = 90) (8.137) 


Differentiating this equation with respect to s,, we have 


dh = ,_—__ tns(4~40) , 
= = 


anise a 8.138 
dsp vee ( ) 


Combining this with equation (8.136) implies that 


dsp 


eae Se 8.139 
=, (8.139) 


which means s, = t by using the initial condition s, = 0 when t = 0. Now we have 


=m.,0 m,(1 — ¢o) 
h(t) = ,/ ———~ erf | —-———-1]. 8.140 
From these solutions, we can calculate the time t* when the porosity decreases to 


the typical transition value of ¢”*, 


_ f= 1-¢ 
— Vin i ey oe 
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and h(t*) = Io, 


Hy = fe ent). (8.142) 


The typical value of h at the transition is 


» ge A OG 
hay ; 14 
Ms 7 (8.143) 
The values of solid velocity u* and effective pressure p* at the transition are re- 
spectively 
; ms(1 ae 0) 
* = h — ——_— 144 
u = (8.144) 
and 
22 1— ¢* 
P= 9ne( 1S - In : 8.145 
p* = rid = bo) f (8.145) 


It is clearly seen that the above solutions are only valid in the top region with a 
depth of IIp. Below this region, the approximation is invalid, and we may expect a 
transition region which joins the regions where ¢ > ¢* and @ < ¢”*. 

In the outer region just above this transition layer as h — z — IIp, we can write 
the solutions approximately as 

on ge +o%(2—h+ Th), 


* 


n= eS PAT), 


= 
aa 


prp'—-(1-¢"*)(e-h+I), (8.146) 


where we have used u, = —p/= and p, = —(1—¢). ¢* is a constant which is now 


determined. Using ¢(h — Ilo, t) = ¢* and mass conservation, we have 
dit oh =0 at z=hA-Tl, (8.147) 


and 
@&=(1-¢")u,-—u"d, at z=h-TL. (8.148) 


Combining these equations, we get 


ozh + (1— $*)u, — u*d. =0, (8.149) 
or 
_ pt—¢) 8 
aia 8) A oot 
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. ¢ / 
from which we can write ¢@* as 


gS, (8.151) 


Transition Layer 
In the transition layer, we define 7 = o = me, 
alg 
m 


thus ¢~ ¢* + ¢" (B+ iL) ~ o* exp[y(B + 7+)|. Therefore, we define 


and put 

2= hh b+]; B<<1, (8.152) 
wv 

~= dp exp[C + map and C=7B, (8.153) 


whence it follows by a matching principle that VU ~ yn as 7 — oo. 


Based on u ~ u* — Z(B + +), we anticipate that u* — ve ~ +. Therefore, we 
put u= v and we have 
W ~ m(u* - a y= Pn =W*- Pon W* = m(u* — = ), (8.154) 
as 1) — 00. 


Using the relations 0, = md,, 0, = 0; — mhdr, the governing equations become 


1 Ww 
—¢*eCtm[(—W, — hW,] + (1— bet )W, — —d*etmU, =0, (8.155) 
m m 
WwW 
— = emp, + (1 — pret) (8.156) 
p= =o, (8.157) 
By choosing C' = —2 Inm to balance the terms in the above equations, we have the 


leading order approximations 


hgW, + (1—¢*)W, =0, (8.158) 
W = -ep,, (8.159) 
p = —=W,,. (8.160) 


These equations are subject to the matching conditions 


* 


D 


Urn, pr p, Wn~wWwWr- lh 8S Oe, (8.161) 
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The last two equations give 
Weve: (8.162) 


Integrating the first equation and using the matching conditions, we have 


* 


how + (1-6 W = (1- 8°)(W* — En) + hd*yn=(1-6W*, (8.163) 
where we have used y = ” a wae oe and u* << 1. Thus we have 
ho* 
=WwW*- U 8.164 
W=W 1_-¢ ( ) 
Combining this equation with (8.162), we finally have 
a l= OW"... 
EU, = (UW - ee ae (8.165) 


This is of the form of a nonlinear oscillator =W,,,, + V’(W) = 0 with the potential 


Viv) =(K-We and K= os (8.166) 


Vets = ie X, (8.167) 


which reaches its maximum at V = K. The only trajectory which can match to a 


solution in 7 <0 is the one with VW — K as n — —oo. Therefore, we define 


1 | ok w* 
We 1 = ee Lea (8.168) 
h¢* 
and we require V — WV, as 7 — —oo. 
Now rewriting equation (8.165) in terms of V and 7, we have 
BW, + V(b) =0, (8.169) 


with a matching condition V, ~ y as 7 — oo. Integrating this equation, we have 


1 1 
Wy? W) = Ey? 8.170 
aun t VE) = 5B’, (8.170) 


which is an energy equation. We also require V, = 0 when V = K, that is 


1 
5a = (8.171) 


CHAPTER 8. PRESSURE SOLUTION CREEP AND VISCOUS COMPACTION 157 


This relation determines V,,. Once we have V,,, we can determine W* from (8.168). 


That is 
_ he*Vy 


We (8.172) 
Combining (8.144) and (8.154), we have 
a gla Oo, 2R a 
k= Tis(5 = my Inm 4 a (8.173) 


Substituting the expressions for p*, W*, we finally have 


hes 1. oh ; (8.174) 
m(1—¢*) 


which determines h and is valid for t > t* or h > Ip. It is clearly seen that the 
leading order approximation of (8.174) is constant with (8.143) when t = ¢*. 
Solution below the transition layer 


From the transition layer, we have W,p — 0, V — Vy as n — —oo. Now going 
back to z rather than 7 and substituting 7 = m[z — h — Hp — B] into the equations 
(8.155)-(8.157), we have 


W 
=p GE, SOS Peer yw, = — et, = 0, 
m 


De (1— Maal 
p= ——EW,. (8.175) 


In order to balance the second and the third equations, we suppose that 


1s 1 1 
m mm 


Ww. (8.176) 


m 
Now the governing equations become approximately 
—¢*U,+ (1— ¢*)W, =0, 
W = —e¥~(1— 9"), 


A 


p=-EW,. (8.177) 
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By using the boundary conditions W, D, V0 6593 f= 1lp— - Inm, we can 
easily write the solutions for these three equations as 


A 


p=V=0, W=-e%9(1-¢*). (8.178) 


The fact that the constant W does not satisfy the boundary condition W =0 at the 
base z = 0 suggests that there should exist a boundary layer near the base although 
its thickness is only of the order of O(+). 

Boundary layer at z = 0 


In the boundary layer at z = 0, V © Uy, then we have approximately 


Pee 2Inm WwW 


ut = (SEY"Eul, —(L— dee), and ba = 6 exp[-—E™ +=], (8.179) 
with a boundary condition 
= Bie 20. (8.180) 
The outer solution is u® = =(f 1 — 0), which implies a far field condition 
Uz 70 as z— oo. The solution is 
s Ox m : (Se yme 
ue = ( D (1 — $00) [e =i) (8.181) 


This completes the solution procedure. 


Comparisons 


The comparison of the solutions in the lower and upper regions (dashed) with numer- 
ical results is shown in Fig. 8.7 for the case of A = 100 and t = 5. Figure 8.8 shows 
the comparison for the basin thickness h(t). A reasonably good agreement verifies 
the validity of the solution procedure. 

We see clearly that in the case of fast compaction (A >> 1), compaction occurs 
throughout the basin, and the basic equilibrium solution near the surface is a near 
parabolic profile of porosity. However, as depth increases, the permeability has de- 
creased sufficiently, and there is a narrower transition region which marks the sharp 


variation of permeability with porosity. More generally, we might therefore expect 


CHAPTER 8. PRESSURE SOLUTION CREEP AND VISCOUS COMPACTION 159 


that in a marine environment where stratigraphic layers cause sudden changes in 
permeability, that clay layers with small permeability may be associated with the 


formation of abnormally high pressures. 


Scaled height: Z 
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Figure 8.7 Comparison of travelling wave solution (8.135), and tran- 
sition solution (8.170) (dashed) with the numerical results (solid) for 
= 100) att=5. 
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Figure 8.8 Comparison of short time solution (8.128) and large time 
solution (8.174) (dashed) with the numerical results (solid) for 4 = 100) 
at t <5. The two dashed curves joins at t = t* (or h = Ig), and there 


is a jump [A]i+ = O(2) which is relatively small. 


m 
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8.5.4 Summary 


For a model of sedimentary basin formation which incorporates viscous compaction 
due to pressure solution, we have been able to derive approximate solutions in the 
distinct limits of slow compaction (A < 1) and fast compaction (A >> 1). When 
\ <1, compaction is limited to a basal boundary layer of thickness O(y/d) (8.106). 
This result is similar to that which occurs for elastic compaction, and is equivalent 
to results obtained in viscous compaction in the asthenosphere (McKenzie 1984). 
When A > 1, compaction occurs throughout the basin, and the basic equilibrium 
solution (which we may call normally pressured, since the pore pressure increases 
hydrostatically) which applies near the surface is a near parabolic profile of porosity 
as a function of depth. 
bg ae (8.182) 
22M, 
this compares with the equilibrium elastic profile, which is exponentially decreasing 
with depth. However, this normally pressured solution is only valid to a depth Ilo, 
given by (8.142), and approximately 
2m, _ (1-o*\)'? 
m= {om (=<) (8.183) 


At this depth, the permeability has decreased sufficiently, and there is a narrower 


transition region which marks the sharp variation of permeability with porosity. No- 
tice also that even if the permeability exponent is not large, so that ¢* is small 
(d* = do exp[—4In J), nevertheless (8.183) implies that the critical depth is still fi- 
nite. Thus the switch from normally pressured to abnormally pressured is predicted 
to occur in any case in a marine environment where stratigraphic layers cause sudden 
changes in permeability and the subsequent formation of abnormally high pressures. 

At greater depths still, cementation begins to occur. As the grain boundaries 
begin to become cemented, pressure solution will decrease, and it can be expected 
that the rheology reverts to an elastic one; from the point of view of the sediments, 
compaction will cease and the medium will become virtually rigid, with pore pressure 
being controlled purely by Darcy flow. Incorporation of these and other processes 


such as diagenesis will form the substance of future developments. 


Chapter 9 


Conclusions 


Conventional studies have treated compaction and diagenesis separately, and the ex- 
perimentally derived nonlinear behaviour of soils is generally studied numerically and 
has never been investigated on a basin scale. One main novelty of the current re- 
search work, based on previous work (chapters 1-2), is to model compaction, thermal 
history (chapters 2-4), unloading (chapter 5), diagenetic (smectite-illite) reactions 
(chapters 6-7) and pressure solution creep (chapter 8) simultaneously in a compact- 
ing 3-D frame. Another novelty is that the nonlinear sediment behaviour including 
phenomena such as hysteresis is treated for the first time on a basin scale, and the 
basinwide response to the unloading from surface erosion is also investigated (chap- 
ter 5). The coupled partial differential equations with a free boundary are solved 
numerically, and analytical solutions are obtained for some geologically interesting 
cases such as rapid and slow sedimentations, diagenesis in the temperature range of 
hydrocarbon generations and basinal response to sediment erosion at the surface in 
a nearly equilibrium state. The well known Athy’s law, derived from field data, can 
be easily obtained from our analytical solutions. 

Based on the pseudo-steady state approximations, the transport-reaction model 
equations of compaction and diagenesis can be simply written in dimensionless form 


as 
a(1 — ¢) 
Ot 


eo =0, (Mass conservation) (9.1) 


Oz 
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ee 
v= — Ak + (1—)]. (Darcy's law) (9.2) 
A constitutive relation is needed to complete this model. In the case of poro-elastic 


compaction, we use an Athy-type relation p = p(@) (chapter 4); while in the case of 


—n Ous 


Oz 


viscous compaction due to pressure solution creep only, we choose p = =o(4) 
(chapter 8). These two different rheological relations result in two quite different 
behaviours of porosity evolution. In the simpler poro-elastic case, we have a single 


non-linear diffusion equation for porosity @ 


Oo Oh ot _ 2,1 Ob 
BE 9, tk o)'| 


@ Oz 
k = ($/¢0)™, m=8, (9.4) 
with a moving boundary described by 
i =1+ dC - 0558 — 1), (9.5) 
and boundary conditions 
i, =o=0 at. 2=0, (9.6) 
b= at z=h, (9.7) 


The analysis in Chapter 4, which was further elaborated by Fowler & Yang (1997), 
showed that the limit \ << 1 (slow compaction) can be simply analysed by means of 
a boundary layer analysis at the sediment base. The more interesting mathematical 
case is when A >> 1 (fast compaction). For sufficiently small times, the porosity 
profile is exponential with depth, corresponding to an equilibrium (long-time) profile. 
However, because of the large exponent m in the permeability law k = (¢/¢9), we 
find that even if A >> 1, the product Ak may become small at sufficiently large 
depths. In this case, the porosity profile consists of an upper part near the surface 
where Ak >> 1 and the equilibrium is attained, and a lower part where Ak << 1, 
and the porosity is higher than equilibrium. Straightforward asymptotic methods are 
difficult to implement because the limit m >> 1 implies exponential asymptotics, 
but we use a hybrid method which appears to correspond accurately to numerical 


computations. 
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For the case of viscous compaction due to pressure solution creep, the equations 


can be simplified as, by taking n = 0, 


a6 8 7 
2g |OONaa. 
= ak | 1 a). 
O 


whose boundary conditions are that 
u=0 on z=0, 


p=0, 6=¢, h=me+u at z=A(t). (9.9) 


The analysis in Chapter 8, which was given in more detail by Fowler & Yang (1998), 
showed that for 41 << 1, compaction is limited to a basal boundary layer of thickness 
O(/2). This result is similar to that which occurs for elastic compaction, and is 
equivalent to results obtained in viscous compaction (McKenzie, 1984; Birchwood & 
Turcotte, 1994). 

When A > 1, compaction occurs throughout the basin, and the basic equilibrium 
solution (which we may call normally pressured, since the pore pressure increases 


hydrostatically) which applies near the surface is a near parabolic profile of porosity 


%, (1 — $0)? 2, 
DO = pag Ve) (9.10) 


this compares with the equilibrium elastic profile, which is exponentially decreasing 
with depth. 
However, this normally pressured solution is only valid to a depth II, given by 


(8.142), and approximately 


_f 2m, (1-¢t\\"” 
nan { i (=) on 


At this depth, the permeability has decreased sufficiently that the hydrostatic balance 


no longer applies, and there is a narrower transition region in which the effective 


pressure drops to near zero, and the porosity profile changes shape. This transition 
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region marks a (relatively sudden) switch from a normally pressured environment to 
one with high pore pressures, and is caused by the sharp variation of permeability 
with porosity. Notice also that even if the permeability exponent is not large, so that 
o* is small (¢* = ¢9 exp[—4In ]), nevertheless (9.11) implies that the critical depth 
is still finite. Thus the switch from normally pressured to abnormally pressured is 
predicted to occur in any case. More generally, we might therefore expect that in a 
marine environment where stratigraphic layers cause sudden changes in permeability, 
that clay layers with small permeability may be associated with the formation of 


abnormally high pressures. 


9.1 Main Conclusions 


A general mathematical model of compaction and diagenesis is presented in this work. 
The coupled non-linear diffusion equations have been solved numerically and several 
asymptotic solutions are given for the cases of geological importance. Asymptotic 


analysis and numerical simulations showed that 


e The processes of diagenesis, temperature and porosity evolution for continuously 
compacting sediments are characterised and controlled by three dimensionless 
parameters A, A, R which relate the sedimentation rate, permeability, heat con- 


ductivity, viscosity, diagenetic reaction rate and heat flux. 


e The present model clearly degenerates to that of Audet & Fowler (1992) by 
setting a = 0, k, = 0 and omitting the temperature equation, or equivalently 


leaving out the parameters A and R by setting them to zero. 


e The parameter \ = 1 defines a transition between slow sedimentation and fast 
sedimentation. Here, the fast and slow are only meaningful relative to the time 
scale for the compaction process. 4 >> 1 corresponds to the case of slow sed- 
imentation or high permeability and \ << 1 corresponds to that of fast sedi- 
mentation or low permeability. The parameter governs the mechanism of the 
excess pressure development of the sedimentary basins. High sedimentation rate 


may cause excess pressures even in basins with moderate permeability. 
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e The effect of variation of sedimentation rates and unloading are investigated 
by solving two sets of equations with a switching condition derived from critical 
state soil mechanics. A very interesting phenomenon arises in the case of constant 
unloading. A downward travelling interface separates the unloading region from 
the loading region. If the system is reloaded, a discontinuous porosity may occur. 
A literature survey suggests that this is the first attempt to investigate unloading 


on a basin scale. 


e The parameter A also defines a transition in a quite similar manner. A << 1 
shows that the temperature solution is dominated by the fast moving bound- 
ary effect of the basin due to fast sedimentation, while A >> 1 shows that the 
sedimentation rate has an negligible influence on the temperature development. 
In the realistic geological environment, it is usually A >> 1, which means that 
the time scale of thermal conduction is much shorter than the time scale of com- 
paction, thus temperature evolution is essentially independent of the compaction 


process as the coupling is very weak. 


e The parameter R, which may be defined in terms of a critical temperature ©,, 
controls the speed of diagenesis and its characteristic effect on compaction. This 
study reveals that mechanical compaction, which is controlled by the strata 
permeability and sedimentation rate, is the most important geological factor in 
porosity reduction and the formation of overpressure. The chemical compaction 
controlled by the diagenesis is of secondary importance in the whole mechanism. 
The first-order dehydration model of diagenesis is a good approximation in the 
sense of only describing the extent of progress of the overall smectite-to-illite 


transformation without much care of its detailed geochemical features. 


e Diagenesis has been successfully modelled as a dissolution-precipitation reaction 
model which can reproduce many essential features of the smectite-to-illite pro- 
cess if the appropriate reaction rate laws are used based on the known physics 


and chemistry from experimental studies. 


e Pressure solution is an effective compaction mechanism as well as a source of 
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cementing material, especially in sandstones and carbonate packstones. Pressure 
solution is successfully modelled as a viscous compaction creep model. Athy’s 


law is replaced by a viscous rheology, and the present model extends earlier work. 


9.2 Future Work 


Based on our current research work, the main further research objectives are 


e Application to practical sedimentary basins. Our new models will be tested 
by real field data, and modifications will be added to the models in order to 
make more realistic modelling. We aim at making reasonable predictions of 
overpressuring before drilling and identification of its precursors by using in situ 
data; and modelling hydrocarbon generation and migration to predict reservoir 
quality. We also aim at applying the present diagenetic reaction model to other 
clay minerals such as quartz production, and geochemical weathering processes 


in the near surface environment. 


e EHxtensions to formulations of new rate laws of natural water-rock systems. One 
unsolved problem in the studies of water-rock interactions is that the laboratory 
data are not directly applicable to field observations. The discrepancies between 
field estimates and laboratory measurements of reaction rates are as large as up 
to four orders of magnitude (Swoboda-Colberg & Drever, 1993; Bitzer, 1996). 
Therefore, we intend to aim at extending our present models to formulate more 


realistic rate laws of water-rock systems in the field. 


e Development to model the nonlinear sediment creep behaviour. Biot theory in 
soil mechanics prescribes a relation between porosity (volume fraction of pore 
space in total volume) and effective pressure (overburden pressure minus pore 
pressure). It is valid for small strains and usually in the one-dimensional case. 
It is not true for large strains. We therefore intend to aim at developing a 
fully nonlinear soil mechanic model to correct the drawbacks of the present Biot 
theory and to reproduce much of the experimentally derived nonlinear features 


of sediments. 
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An eventual aim of the model development will be the production of a code which 
can solve the compaction problem which includes most of the known physics and 


chemistry, in a three-dimensional environment. 
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